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FOREWORD 


This research was conducted by The Pennsylvania State University for 
the National Aeronautics and Space Administration, Georgy C. Marshall 
Space FI i^t Center, Huntsville, Alabama. Dr. George H. Fichtl of the 
Aerospace Environment Division, Space Sciences Laboratory, was the 
scientific monitor. Professor John A. Dutton, Department of Meteorology, 
The Pennsylvania State University, was the principal investigator. The 
support for this research was provided by Mr. John Etiders of the Aviation 
Safety Technology Branch, Office of Aeronautics and Space Technology, 

NASA Headquarters. 

The research reported herein was motivated by the need in the aviation 
community for advanced turbulence simulation schemes that provide 
simulated turbulence time histories for flight simulation applications 
which satisfy not only the known second-order statistics of the atmos- 
pheric turbulence as in the state-of-the-art turbulence simulation 
schemes now being used by the aviation community, but in addition 
include information concerning higher-order statistics as manifest by 
the non-Gaussian nature of atmospheric turbulence. In this context, 
it is believed that the work presented herein represents a major 
step forward in stochastic process theory as related to atmospheric 
turbulence simulation. However, it should be noted that the results 
of this report, when used alone, are only applicable to the case in 
which an aeronautical system is completely immersed in the gusts; 
i.e., absence of gust gradient effects. In this regard, additional 
work is under way to develop gust gradient simulation schemes whereby 
the gust velocity time histories generated by the technique reported 
herein and the gust gradients simulated by the techniques now under 
development are applied with a Taylor series expansion (truncated 
at the first-order term) of the gust velocity vector field about the 
vehicle center of gravity, so that the gust environment at any 
point on the vehicle can be generated during a flight simulation. 
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SUMMARY 


The method of simulating a turbulent time series by filtering a 
white noise series Is consolidated and extended. The development of 
linear filters from empirical spectra Is expanded for forms based on 
boundary layer similarity under stratification and for generalized 
spectral shapes. Some properties of the filters under various stratifi- 
cations, heights and viscosities are examined. The method of linear 
simulation is extended to multi- component processes by the dlagonal- 
Izatlon of the spectral matrix, spectral factorization of the 
eigenvalues, followed by a rotation Involving a special unitary trans- 
formation. Results Indicate that the addition of a cross-response 
increases the total response. 

The method of simulation by filtering is extended to several non- 
linear, non-Gausslan models. These models are based on ad hoc approxi- 
mations of the kernel Interactions. It Is found that the method of 
separable kernels for a representation of velocity is inappropriate 
for simulating the characteristic Inertial transfer of turbulent energy. 
However, the separable kernel representation of acceleration better 
approximates the energy transfer in the viscous subrange for suf- 
ficiently small Reynolds number. 

An evaluation of the linear and non-linear models, with computa- 
tions carried out In phase space. Is Included. Because the non-linear 
simulation method requires a more precisely Gaussian stimulating process 
than Is commonly available, special generative techniques were developed 
and examined. Results Indicate that non-linear simulations will re- 
quire large arrays of very nearly white, Gaussian noise In applications. 
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1.0 INTRODUCTIOK 


The Increasing sophistication of design techniques In a variety of 
engineering and environmental applications requires the simulation of 
the statistical structure of wind gusts. As ioan~made structures have 
been made larger and more flexible. It has become Imperative that the 
effects of the wind, both as a static and fluctuating force, be In- 
corporated In structural and economic designs. In addition. In 
response to Increasing concern for the quality of the environment, the 
need to understand the wlnd~drlven diffusion mechanism distributing 
alr-bome pollutants through the atmospheric environment has In- 
creased. Fortunately, there exists a wealth of Information about the 
statistical structure of the wind, particularly near the ground. How- 
ever, an Important difficulty lies In Incorporating such Information 
Into applications In a manner that Is at once practical and realistic 
as It conveys Important aspects of the meteorological dynamics of the 
problem. 

1.1 Statement of the Problem 

In this study, the method of Incorporating the statistical 
structure of the turbulent wind field near the ground Into applications 
Is consolidated and expanded. 

The linear spectral representation of turbulence has provided a 
tiseful Interface between the meteorologist and the engineer. Conse- 
quently, Its properties, successes and failures are well known. The 
Interest of the meteorologist has been generally centered on providing 
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the best possible estimates about the structure of the wind field, both 
as It relates to the vertical distribution of the averaged wind and to 
the spectral distribution of the fluctuating wind. 

Models to simulate turbulence which require the specification of 
the mean wind with height and thermal stratification vary In their 
sophistication and reliance on the principles of boundary layer simi- 
larity. As a result, the description of the vertical profile of wind 
has tended to be a potpourri of empirical relationships and approxima- 
tions. Because the numerous parameters which characterize turbulence 
In the atmospheric boundary layer have often not been measured 
simultaneously, the applicability of some empirical results Is unknown. 

Accordingly, there Is a need to consolidate aspects of the verti- 
cal structure of the wind, in order to have them consistent with known 
similarity properties of the flow, and to be able to Incorporate fur- 
ther results as they become available. 

Another crucial feature of simulation models of turbulence near 
the ground Is the approximation used for the spectral distribution of 
the variance In the fluctuating components of the wind, particularly, 

In the range of scales of size equal to or less than the distance, z, 
from the ground. Simulation models which lead to modeled realizations 

of the turbulence have required the so-called Dryden spectral form, 

-2 

which for sufficiently small scales varies as k , where k la the wave 

number. This spectral form has no basis in theory or observations and 

has been chosen only for Its analytical properties. Other applications 

not requiring modeled realizations have been based on the von Karman 

- 5/3 

spectral form which tends to k , for kz » 1, In accordance with 
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the well-known properties of the Kolmogorov Inertial subrange. However, 
the von Karmw spectrum requires the specification of a length with 
which to characterize the bandwidth of the spectrum, but this length 
parameter bears no known dynamical relationship to the structural 
properties of the turbulence in the atmospheric boundary layer. As a 
result, the implied spectral dynamics of von Kdrman's model cannot be 
deteinnined by recourse to theoretical considerations of the boundary 
layer. The determination of such lengths must of necessity be made 
empirically. It remains to determine alternative spectral forms which 
are dynamically consistent with, say, the vertical profile of wind 
and dissipation determined by similarity arguments. 

It is also well known that the concept of a linear representation 
of turbulence in terms of a Gaussian, white noise process is inconsis- 
tent with the observed non-Gaussian and non-linear structure of the 
turbulence. In particular, linear Gaussian models are inadequate for 
the simulation of the large gust structure. Therefore, an extension 
of the representation of turbulence is considered in this study which 
systematically incorporates some basic properties of the non-linear 
and non-Gaussian probabilistic structure of the turbulence. The 
mathematical formulation of this extension is most conveniently based 
on a functional series expansion in terms of the simple and convenient 
Gaussian, white noise process — the same as is used in linear modeling. 

The method of functional representation will be shown to lead 
naturally to a concept of a discrete gust form. As such, the method 
of the representation is superior to other discrete gust models where 
mathematically convenient, ad hoc forms are specified. Because here 
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the discrete gust form Is a derived property of the process, a consist 
tent formulation for the statistical structure of the turbulence In 
the boundary layer allows for a systematic analysis of this discrete 
gust form for various thermal stratifications and Reynolds number ef- 
fects. 

In summary, It Is the purpose of this study to consolidate the 
simulation modeling of turbulence In the boundary layer In terms of 
boundary layer similarity principles and empirical results, and to 
extend the modeling for some aspects of the non-llnear and non-Gausslan 
structure of the turbulence there. It Is also the Intention of this 
work to Identify some properties of the discrete gust form structure 
of the modeled turbulence. 

1.2 Importance of the Problem 

The study of the structural effects of turbulence In the earth^s 
boundary layer divides naturally Into three main streams of research 
and development (1) the collection and assimilation of turbulence data; 
(2) the theoretical modeling of the statistical and dynamical nature 
of the turbulence; and (3) the developm^t of methods by which to 
specify the response, structural or otherwise, to the turbulence. 

While a considerable effort has been made to measure, describe and 
model the flow field throughout the planetary boundary layer, the 
development of methods of application of this accumulated wealth of 
Information has not been rapid. 

The need for applied models of turbulence near the ground Is 
ubiquitous. The need Is perhaps most severe in the design and opera- 
tion of aircraft and other aerospace vehicles. In aeronautical 
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operations, cross-winds and wind shear nay on occasion present a hazard 
In the approach and landing stages of a flights In addition to ques- 
tions of structural Integrity and passenger comfort, the principal 
concert of the designer Is that the pilot may lose control as the air- 
craft Is accelerated, or may aggravate the situation by Initiating the 
wrong corrective procedure. The control problem Is compounded for 
VSTOL aircraft and helicopters whose lift characteristics are more 
sensitive to the direction of air flow relative to the lifting surfaces. 
More recently, the need has developed to study the response of rockets 
to turbulence during launch because the wind fluctuations affect the 
stability and navigation of the vehicle. 

The concern about wind effects is shared by many other engineers. 
Turbulent buffeting of surface structures, particularly of large 
flexible bridges and office towers, must be considered at the design 
stage. Further, the design of surface transportation systems and 
vehicles also requires a specification of the range of probabilities 
of significant wind events and a method of estimating the response or 
result. Of Interest, particularly In large urban areas. Is the effect 
of turbulence In the dispersion of pollutants in the atmosphere. Yet 
another Important area requiring the modeling of the wind field Is 
that of water wave formation and maintenance. In as much as the wave 
environment effects over-water transportation, recreation, and the 
dispersion of pollutants In or on the water. 

1.3 General Characteristics of Simulation Models 

Of course, not all requirements for simulating the statistical 
structure of the wind field can be met by developing a particular model. 
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For example, the methods used to simulate the dynamics of turbulence 
(Kralchnan, 1965; Herring, 1966; Deardorff, 1972a), while useful In 
testing the consequences of the approximations characterizing each 
model, are either not theoretically compatible with the Inhomogeneous 
structure of turbulence In the atmospheric boundary layer or else are 
inq>ractlcal to Implement. 

The method found most suitable by engineers for simulating 
turbulence Is the so-called spectral filtering, or forcing, technique 
pioneered by Llepmann (1954) . The spectral filtering model essentially 
characterizes the response of an aircraft or structure, or any process 
driven by turbulence, as a signal derived by filtering a sequence of 
pulses uncorrelated sequentially whose amplitudes are derived 
probabilistically from a Gaussian distribution* The desirable property 
of the latter process is the constant spectrum, which, from optics. Is 
referred to as a white spectrum. Because of the wide variety of ap- 
llcatlons (Houbolt, 1973), the method remains a useful technique. Its 
success to a considerable degree Is attributable to Its simplicity. 

The spectral filtering method, as the name implies, Is based on the 
use of the spectrum of the atmospheric turbulence to characterize the 
flow field. The characterization of the spectrum over scales Important 
to the application In turn requires specifying the spectral form as 
well as Its variance and its bandwidth. Variations exist In the 
representation of the spectriun and its controlling variables (Teunlssen, 
1970), and are discussed In more detail later. The mathematical 
details of the method are also postponed for later consideration. 

The simulated turbulence resulting from the linear filtering of a 
Gaussian process Is Itself Gaussian. However, turbulence Is not a 
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Gaussian process (Dutton, 1970) . In order to produce a non’-Gausalan 
process which simulates turbulence from a white, Gaussian process 
requires non-*llnear filters. The Introduction of non-llnear filters to 
synthesize turbulence Is a recent development (Reeves, 1969; Kurkowskl, 
et al., 1971; Gerlach, et al., 1973). However, the non-linear models 
which have been developed suffer from a lack of consistency in terms 
of the observed, non-llnear properties of boundary layer turbulence. 
Therefore, what Is needed Is a systematic methodology by which to 
Introduce the observations of the non-llnearlty of atmospheric turbu- 
lence Into the filtering method. 

1.4 Structure of the Atmospheric Boundary Layer 

For a steady, horizontally homogeneous mean flow In the boundary 
layer, sufficiently near the ground, the vertical variation of the 
turbulence fluxes Is negligible (Blackadar and Tennekes, 1968); 

In particular 

t(z)/p =» - uw = u.^ (1.4.1.) 

o ** 


and 


H(z)/p^ Cp " w0 - - u^ T^ (1.4.2) 

are independent of height. In (1.4.1) and (1.4.2), x/p^ and H(z)/p^ Cp 
are the specific momentum and heat fluxes respectively; p^ Is the 
density of the air and c the specific heat at constant pressure. 
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Equations (1.4.1) and (1.4.2) also serve as definitions for the 
characteristic velocity, and temperature T^. According to the 
hypothesis of Mbnln and Oboukhov (1954) the structure of the mean wind 
shear and temperature gradient (sufficiently near the ground so that 
Inertial effects due to the earth's rotation are unimportant) can be 
derived on the basis of dimensional arguments. That Is, the mean 
shear and temperature gradient are given by 


and 


dz 




(1.A.4) 


where L, the Monln*-Oboukhov length. Is 


L 


u . T 
* o 

g K T 


* 


(1.4.5) 


In (1.4.3) to (1.4.5), U Is the mean wind speed, 6 the mean 
potential temperature, the depth averaged boundary temperature, ic 
(a similarity parameter) is von Karman's constant, and and (j)^ are 
the similarity functions for the shear and the terperature gradient. 
It Is convenient In what follows to define z/L by 


C " z/L 


(1.4.6) 
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In fact, the Mbnln-Oboukhov hypothesla states that all statistics of 

the turbulence sufficiently near the ground. In dlabatlc situations, 

for Ideal steady, horizontally homogeneous flows, become functions of 

C only. If velocities are scaled with u^, temperatures by and 

lengths by kz. Accordingly the moments 1^, of the probabilistic density 
lltl 

function of the 1 turbulence velocity components are also functions 
of C* 



(1.4.7) 


The spectral distribution, of variance or covariance over the range 

of wave numbers, k, for which there Is any shear or buoyantly Induced 
turbulence, becomes, under the appropriate scaling. 



Gy (f.C) 


(1.4.8) 


In (1.4.8), represents the appropriate scaling variable, u^ or T^, 
and f, where 


f * KZ k 


(1.4.9) 


represents a normalized wave number. The use of subscripts 1 through 
3 assumes the standard meteorological usage (Lumley and Fanofsky, 

1964) and the subscript 6 refers to temperature. 

For scales, f » 1, such that the turbulence becomes asymptotically 
Independent of the details of the mechanisms generating the turbulence 
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(Tennekes and Lumley, 1972, Chapter 8), the spectra, ^11 are similar. 
In the sense of dimensional analysis, under the hypothesis of 
Kolmogorov (1941), so that 

Gji(£) - (1.4.10) 

Specifically, for a sufficiently large Reynolds number, defined by 

Rey Kz u^/v (1.4.11) 

where V Is the kinematic viscosity, 

♦ll ■ *22 ■ *33 ■ (1.4.12) 


and 


♦ee - 


-1/3 


a) 


(1.4.13) 


The functions <)> and ({) represent the similarity functions under 
X ^ 

Monln~0boukhov scaling, for the dissipation rates of kinetic energy, 
€, and temperature variance, x» defined by 


<l>p(C) • Kze/u^' 


(1.4.14) 


and 


4>^(0 - Kzx/u*T*^ 


(1.4.15) 
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A considerable effort has been made by many Investigators to 

Identify the similarity structure, both of the low order moments 

of turbulence near the ground, under diabetic conditions, and of the 

spectral and co-spectral forms. An excellent summary and review Is 

provided by several authors In a monograph edited by Haugen (1973) , 

and will not be duplicated here. A summary of empirical forms for 

♦ * as functions of C* are given In Appendix A, 

tt n E X ^ 

as well as empirical spectra and co-spectra, 6^^ , as functions of 
( and f. 

In order to utilize the similarity relationships of Monln and 
Oboukhov, It Is necessary to estimate the factors u^ and by an 
Independent method. For a steady, horizontally homogeneous boundary 
layer, Kazanskl and Monln (1961) derived the resistance laws 


ln(G/fz ) - 8 + In (G/u^) + (k^ 


(1.A.16) 


and 


sin “ ■ “ 7^ (1.4.17) 

In (1.4.16) and (1.4.17), G Is the geos trophic wind modulus and f Is 
the Coriolis parameter (* 2(2 sin where (2 Is the earth's angular 
velocity and Is latitude) , a Is the angle between the direction of 
the surface stress and the geos trophic velocity, and z^ Is the surface 
roughness. The functions, A and 8, are similarity functions, which 
for dlabatlc conditions are hypothesized to be universal In the 
stability parameter, y given by 
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U - h/L (1.4.18) 

where h Is the so-called Ekman height 

h - K \ijt (1.4.19) 

characteristic of the depth of the boundary layer under neutral 
stratification. An analogous development by Zllltlnkevlch and 
Challkov (1968) for the transfer of heat across a turbulent boundary 
layer Is given by 

Ae/T^ « P^[ln(u*/f z^) - C(y)] (1.4.20) 

where A6 Is defined as the potential temperature difference between 
the surface and the level where the flow is geos trophic* is the 
turbulent Prandtl number (- 0.7) under near-neutral conditions and C 
Is a universal function of y. From empirical formulations for A(y)* 
8(y) and C(y), (Kerman, 1974a), It is possible to construct algorithms 
for the momentum drag coefficient, u^/G and the 'thermal drag' 
coefficient, T^/A6, as functions of the dimensionless parameters, Ro 
and S , given by 

s 

Ro - G/fz (1.4.21) 

o 


and 
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S - g/T (A0/G£) (1.4.22) 

o 

Accordingly* and can be derived Immediately given the ex- 
ternalf controlling variables of the problem — G, and A0. 

The combined similarity theories of Monln and Oboukhov and the 
resistance laws make It possible to estimate the turbulent statistics 
at a given height, z. In the constant flux layer, at a particular loca- 
tion and time, having estimated G, z^ and A6. 

It Is emphasized here that since the basic meteorological dynamics, 
as conveyed by the similarity theories, are self-consistent with the 
empirical representations (such as those given In Appendix A) , It Is 
pointless to Introduce additional variables through ^ hoc models of 
the spectra. That Is, In models such as von Karman^s (Telnlssen, 1970, 
p. 40) the scaling lengths there are not linked dynamically to 
similarity theories. In fact. It Is often observed that a form other 
than the von Karman spectral form, may be appropriate for representing 
the large scale structure of the spectrum. For these reasons, only 
spectra based on direct observations are considered In the models 
developed here. 
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2.0 FORMULATION OF THE SIMULATION MODEL 

In order to extend the methods of simulating turbulence to 
processes with rather general spectral and non-Gausslan characteristics. 
It Is necessary to develop a suitable mathematical structure, explore 
some of Its properties, and determine whether Its application Is 
practical. This chapter outlines a particular functional representa- 
tion for turbulence that allows for a systematic development based on a 
Gaussian process. We also consider here the method of spectral 
factorization for calculating a filter for those cases In which a linear 
sub-process can be Identified. 

The concept of a functional representation (Wiener, 1958) of 
turbulence arises from a picture of a turbulent velocity field as the 
result of random Impulses. The process of generating a response to a 
stimulus Is equivalent to a black-box process. The triplet of Input, 
black-box. and response are the characteristic elements of a mathe- 
matical identify called a system. It is the intent of this research to 
consolidate and extend present system representations which use random 
pulses to produce a response resembling turbulence. 

2.1 Functional Representation of Turbulence 

The mathematical formulation used here for analysis of a system is 
based on the original work of Volterra (1930) . A functional transforma- 
tion Is defined as the operation which transforms a function to a 
number — an operation such as a definite integral. In the application 
of functionals to systems, the black-box Is modeled by the functional 
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transformation which maps the temporal function representing the evolved 
history of the Input to the present (single) value of the output. If 
an Input function, C(t), as evolved up to time t. Is represented as a 
point In a space of functions, then the functional transformation maps 
the Input Into a point In a new space of output functions, say, y(t). 

As ^(t) varies with the parameter, t, then the transformation will map 
^(t) to a varying output, y(t). Volterra showed that a continuous 
functional (transformation) could be uniquely approximated by a poly- 
nomial series of functionals given by 


y(t) 


K" 



K^(t;t^) 5(t^)dt^ 


+ 


rr 

—00 —00 


K^(t;tj^,t2) ^(t^) dt^ dt2 

( 2 . 1 . 1 ) 

iOO mQO fOO 

+ 111 5(tj^) ^(tj) ^(tj) 

—00 —00 —00 


dtj_ dtj dt3 + . . . 


(All Integrals hereafter will have a range (-<», <») unless otherwise 
specified.) The functions, K^, are referred to as the kernels of the 
representation. A physically realizable situation. In which the trans- 
formation of any signal can act only on the past of an Input, requires 
that 


K (t { tj^ , t2 , . . . t^) “ 0 ^^1 ^ 


( 2 . 1 . 2 ) 
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For a horlzonally homogeneous turbulent flow whose statistics are ad~ 
vected according to Taylor's hypothesis » the turbulent process Is also 
ten^orally Invariant, or stationary, so that 

K^(t;tj^,t2,...t^) - K^(t-tj^,t-t2 ,t-t^) (2.1.3) 

In addition to determining the transformation of the Input to the 
output functions, the kernels also determine the statistics of the 
output from those of the Input. A simple, linear, temporally Invariant 
system Is represented by 

y(t) - j K^(t-x) C(T) dT (2.1.4) 

For convenience, we take the Input to have zero mean (x ■ 0). 

Formation of statistical averages of the Input and output leads to 

y(t) y(t+T) - Ryy(f)"Jj dXj^ dXj (2.1.5) 

and 

y(t) 5(t+X) - Ryg(T) - J K^(Tj) Rg^(X-Tj) dTj^ , (2.1.6) 

In (2.1.5), the output variance, represented as a transforma- 

tion of the Input variance, hv The utility of using a whlte-nolse 
Input process, defined by 


R^^(t) - 6(t) 


(2.1.7) 
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(6 is the Dirac delta function) Is shown by substituting (2.1.7) In 
(2.1.5) and (2.1.6) to form 

Ryy(T) - I dTj^ (2.1.8) 

R g(T) - K^(T) (2.1.9) 

1 ! 

The kernel, K , Is therefore derivable by cross-correlating the 
Input and output, and the kernels for the functional transformation 
of the white noise process share variance properties with the output. 

It Is lnq>ortant to note that while the above technique of 
determining the kernel of a linear system Is used In many fields of 
engineering, here the supposition that both the Input pulse and output 
response are available for correlation is not valid. The Input forcing 
mechanism, as represented by the functional transformation (2.1.1), Is 
Internal to the fluid and not measurable. This makes the problem of 
determining the kernels In practice more complicated than the usual 
situation where both input and output are available simultaneously. 

The significance of (2.1.8) Is better seen In a spectral 
representation. Consider a Fourier transform defined by 

f(t) • ^ j f(co)e“^*^ du) (2.1.10) 

applied to (2.1.4). The result Is given by 

y(a>) - K^(ai) C(w) 


( 2 . 1 . 11 ) 
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which for the white-noise property 

lil"- 1 


( 2 . 1 . 12 ) 


becomes 


^ 0 '^1 9 

|y| ■ (2.1.13) 

Thus the spectrum of the kernel is equal to that of the output for a 
linear system. Equation (2.1.13) forms the basis for many of the 
applications of the linear simulation of turbulence. A major question 
remaining, then, is how to determine K^, given that its spectrum is 
that of the turbulence. A recent computational development is 
explored in Section 2.4 to determine for a relatively general class 
of spectra. 

An important consideration in developing a simulation model is the 
ease of application. Farente (1970) has outlined the method of treating 
Interacting systems. In most applications of turbulence models, the 
simulated turbulence is used in turn to stimulate a system representing 
a structure or perhaps another geophysical process. By the algebra 
of functionals (Farente) , the final response statistics are derivable 
from those of the turbulence without recourse to actually generating 
simulated realizations. Of course, such a consideration is basic to 
linear filtering, but it is useful to note that its application is also 
valid with non-linear simulations such as discussed in Chapter 5.0. 

As shown by Barrett (1963), the functional expansion of (2.1.1) 
can be made more efficient by an orthogonalization of the basis, or 
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input process, 'Hie question is what process to use for maximum 

efficiency. The philosophy of this study is that in cases where the 
turbulent process may be considered nearly Gaussian, the obvious basis 
to use to expand the flow field is one based on orthogonallzation 
relative to the Gaussian process — the Hermlte polynomials. Some 
questions about the conditions under which such an expression is likely 
to be successful are discussed next. 

2.2 Probabilistic Structure of Surface Turbulence 

For many years (Batchelor, 1953) the one-point turbulent velocity 
probability density function (p.d.f.) was observed to be indistinguish- 
able from a Gaussian distribution. Stewart (1951) was the first to 
establish the pronounced non-Gauss Ian structure of turbulence with de- 
creasing scale. Further investigations of the moment distributions 
over scale (Frenkiel and Klebanoff, 1967) confirm the converse of 
Stewart's work — that there is a quasi-Gauss ian structure at scales 
commensurate with the energy containing sizes. 

Arguments concerning Gaussian structure are not extendable to the 
joint p.d.f. of two velocities at neighboring points because of the non* 
linear effects (Batchelor, Chapter VIII) which led to an Inertial 
transfer of energy across wavenumbers. Because non-linear interactions 
within the turbulence Increase with decreasing scale up to the viscous 
limit, a resulting increase of non-Gausslan characteristics with de- 
creasing scale is to be expected. The probability distribution of the 
dissipation rate 


e - 15V (|i)2 


( 2 . 2 . 1 ) 
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for Isotropic turbulence (or equivalently local accelerations) Is a 
convenient measure of the non~llnear (and non-Gausslan) structure over 
a vide wavenumber region of engineering concern. According to Kolmogorov 
(1962) , Oboulchov (1962) and Grant et al. (1962) , the equilibrium struc~ 
ture Implied In Kolmogorov's original similarity result required a 
refinement to a more local, volume-averaged dissipation rate, < e >. 

It has been suggested (Gurvlch and Yaglom, 1967) that < e > has a log- 
normal distribution. This prediction has been disputed by Tennekes and 
Wyngaard (1972) and Gibson and Masiello (1971) on the basis of experi- 
mental data taken at a very large Reynolds number. At present, the 
only workable hypotheses on the probability structure appear to be 
empirical (Tennekes and Wyngaard, 1972; Frenklel and Kelbanoff, 1967). 

In siumnary, it Is reasonable to attempt a simulation of the energy 
containing structure of surface layer turbulence in terms of a quasl- 
Gausslan process. The fact that the observed structure of the surface 
layer turbulent velocity field Is nearly Gaussian (Appendix A), as 
expected by the preceding discussion. Is encouraging for modeling 
purposes. However, from the discussions of the strongly non-linear 
spectral region. It Is concluded that an expansion about a Gaussian 
process at scales much smaller than the energy generative region Is 
limited. 

2.3 Wlener-Hermlte Functional Representation 

The orthogonal functional Hermlte polynomials based on Input 
realizations, ^(t), drawn from a white, Gaussian, stationary process 
are (Barrett, 1963) 
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H„(0 - 1 (2.3.1) 

o 

- e(tj^) (2.3.2) 

H2(tj^.t2;C) - 5(tj) e(tj) - ^(tj-tj) (2.3.3) 

“3^*1'*2’^3*^^ " ^^*^2^ ^^^3^ " ^^‘l^ S(t2"'=3) 

(2.3.4) 

- 5(tj) «(t3"t^) - C(t^) ^(ti-tj) 


The Wlener-Hemlte (hereafter referred to W*>H) representation of a 
velocity component, u(t) , by a white, Gaussian, stationary process Is 
given by 


u(t) - j K^(t-t^) H^(t^) dt^ 

+ II K^(t-t^,t-t2) H2<tj^,t2) dt^ dt2 (2.3.5) 

+ III K^(t-t^,t-t2*t-tj) H2(tj^,t2,t2> dt^ dt^ dt^ + ... 


where both Input, ^(t), and output, u(t), are understood to have a mean 
of zero. An equivalent representation follows from a Fourier 
transformation (2.1.10) of (2.3.5) 

u(o)) ■ a (oj) Hj^(o)) + I K^(a)j^,ai-o)^) H 2 (u)j^,i»>-Wj^) du)^ 

(2.3.6) 

+ 


• • • 
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The orthogonality conditions for the Hermlte-polynomlals given by 


H^(t 2 ;C) - 1 




H2<t3,t^;€) = 6(t3~t3) 6(t2-t^) 

+ 6(t^-t^) 6(t2-t3) 


(2.3.7) 

(2.3.8) 


(2.3.9) 


or, their Fourier transformed equivalent. 


•X A 

H^(<u^) “ 6 (o)^‘h 02 ) 

"A" ■ A 

H2(o)^,W2) H2 (0)3,0)^) - 6(0)3+103) 6(0)2+10^) 

+ 6(0)3+10^) 6(o)2+t03) 


(2.3.10) 


(2.3.11) 


considerably reduce the con^lexlty of computing the output statistics. 
For example, by utilizing the orthogonality conditions, the expression 
for the spectral density of the u-process, defined by 

A ' ■ A 

u(o)3) 0(0)2) " ^u^'^l^ 6(0)3+102) (2.3.12) 


Is given by 
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- |k^( 0))|^ +i j |k^ + ... (2.3.13) 

Equation (2.3*13) expresses the decomposition of the spectrum of the 
process Into a sequence of positive definite contributions. The 
positive definiteness of the W-H representation Is a desirable feature 
of the method. Not every moment expansion scheme (Ogura, 1963; 
Deardorff, 1972b) can guarantee such a property, and in this respect 
these other methods contain basic inconsistencies. The implicit 
assumption made In using the W-H representation as a representation of 
nearly Gaussian process Is that the contributions to the spectrum will 
tend to concentrate the variance in the low order terms of (2.3.13). 
Because the Hermite polynomials, and hence individual Hermlte 
functionals of the expansion (2.3.5) and (2.3.6) are orthogonal, the 
truncation 


» 0 1 > 2 


(2.3.14) 


Is consistent with the well-known result for linear white, Gaussian 
forcing 


(})(a)) « |K^(w)|^ 


(2.3.15) 


Successive moment expansions arising from (2.3,5) or (2.3.6) become 


progressively more complicated. For example, the skewness (or bi- 
spectrum) — truncated to second order terms is given by 
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u(o)j^) u(o)2) U (a)^-h02) • 47T K'*'(a)j^) {K'*’(a)2) K (u)j^,aj2) 

/\ A ^ A 

+ K^(ujj^) (o)j^-H02) (-a)^,o)j^-Ki)2) 

+ K^(-a)j^,0)^-H02) (2.3.16) 

1 ^2 ^2 ^2 

+ - [K (w^-*p,p) K (a)^+a)2“P*P-WjL) K (p,w^-Hi)2-p) 

+ K (a32“P»p) K (ajj^-hD2“P,p-W2) K (p,Wj^’H02-p) ]}dp + ... 

Equation (2.3.16) Indicates several other features of the W~-H 
representation. In general, a description of turbulence with an 
infinity of moments is equivalent to a description with an infinity 
of kernels. Also higher moments representing the non-Gauss ian structure 
are characterized by interactions between the K^, or equivalently among 
a hierarchy of non-linearities. It is possible, at least in principle, 
to recover one set of statistics from another by solving the (infinite) 
set of coupled integral equations. 

Some simplification is obviously needed. The truncation of the W-H 
expansion is yet another case in which the closure problem of turbulence 
must be faced. Attempts to determine the kernels dynamically (Meecham 
and Siegel, 1964) have been shown to be inappropriate (Orszag and 
Bissonnette, 1967). Attempts to produce the equivalent of a stimulation 
technique (George, 1959; Dutton, 1970) and correlate the input and output 
are not applicable. The method of Robinson (1967a, b) based on Wiener's 
original work as a method of determining the kernels of a linear (or 
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equivalent linear) system, by predictive decomposition Is unwleldly 
and time consuming. 

For a W-H representation of a nearly Gaussian process, where the 
non‘>llnearltles are weak. It Is advantageous to represent higher order 
kernels In terms of lower order kernels (Chapter 5.0). The result Is 
an expansion about linearity, and the determination of the kernels 
reduces to calculating just The next section outlines a method of 

solving for for a generalized class of spectra which Is compact 
and computationally efficient. 


2.4 Spectral Factorization 

As a demonstration of the baalc features of spectral factorization, 
consider the first order linear system driven by white noise defined by 
the differential equation 

y(t) - 5(t) (2.4.1) 


where y(t) Is the response to the white noise ^(t), and T characterizes 
the response time. For convenience, let us scale the problem so that 
T ■■ 1. This equation is often employed (for example, Skelton, 1968) to 
describe aircraft response to turbulence. The solution of (2.4.1.) Is 
given by 


y(t) 


e ^(t-tj^) dt^ 

0 


(2.4.2) 
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The kernel, K^, is given by e for t > 0, The Fourier transform of 
(2.4.1) is 

y(iw-l) - i (2.4.3) 

or, from (2.1.11) 

K^(o 3) - (io)-l)”^ (2.4.4) 

The spectrum follows frpm (2.4.4) as 
4>y(w) - y y - K K-^ 

- (ii0+l)“^ (-ioH-l)’'^ (2.4.5) 

- (ojVl)”^ 

The spectral factorization problem is the inverse problem. Given 

the spectrum, (p * and the fact it was derived from a white, Gaussian 

'^l 1 

stationary process, find K and K . In the above example the spectral 

factors, 1+lU) and l-lo), are well-known and derivable analytically in 

several ways. In fact, applications with this spectrum (often referred 

to as a Dryden spectrum) have been made simply to utilize the known 

spectral factors and simple form of the kernel (even though it is known 

-5/3 

that the turbulence spectrum varies as o) . Such an assumption is 


useful because the analytical determination of the factors of other 
spectra, such as von Karmim's spectrum, is prohibitive. 
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However finding K from Is not a unique process, as the problem 
Is stated above. There are an Infinity of functions, y(t), which one 
could study, either by observation or simulation which would have the 
property 


y y - K K - <|>y (2.4.6) 


''I 

which la the only defining property of K given. It Is necessary 

A 

therefore to distinguish KT from any y which has the same spectrum. 

'^1 1 

As a demonstration of the defining characteristics of K and K 


consider a simpler problem where the Input is a single pulse. 


C(t) - 6(t) 


(2.4.7) 


From (2.4.2) the output response Is 
y(t) - K^(t) (t > 0) 


(2.4.8) 


that Is, the kernel function Is the response for a single pulse. (Ck>n> 
sequently, Is referred to. In what follows, as a kernel or a response 
function.) The function, K^(t), of (2.4.8) Is, from (2.4,2), 


K^(t) - e~‘ t > 0 


(2.4.9) 


0 


t < 0 
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so that for a pulse (S(t)) at t » 0, the response Is Instantaneously 
at Its maxlmun after the Impulse begins. This property Is referred to 
as the minimum phase or delay characteristic, because the modeled physl** 
cal system responds with the minimum possible delay to a change In the 
Input. Physically, of course, a finite delay would be required before 
the system achieved Its maximum response. It Is noted that the modeled 
process (2.4.1) Is free from frictional effects which Intuitively we 
would expect to delay the Initial response. 

The spectral factorization process using the minimum delay criterion 
was also shown by Bode and Shannon (1950) to be equivalent to determining 
a function with a given modulus (spectrum) with zeroes confined to one 
half of the complex plane. The general factorization problem was 
solved by Kolmogorov and Is discussed In detail by Doob (1933) . 

In what follows, the development is heuristic. Also, because the 
remaining development and computations will necessarily be In discrete, 
tabulated form, the formulation Is given In an equivalent, discrete 

/s 

representation. That Is, K (dropping the super-script 1 for con- 
venience) Is re-deflned as 

A 00 

K((d) - E k (2.<k.l0) 

t-0 ‘ 

where Is the kernel tabulated over the Index t and the summation 
limits reflect the condition that 



1 < 0 


(2.4.11) 
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The phase characteristic, H(o)), of the Fourier transform of the kernel. 
Is defined by 


* K(10) 



(u) e 


10(0)) 


(2.4.12) 


Consider a Fourier transformation of the logarithm of K, given by 


In K - Z L e 
t-0 ^ 


lo)t 


(2.4.13) 


L + Z L cos o)t + 1 Z sin cat 


t-1 


t«l 


The one-'Slded nature of (2.4.13) results from the equivalence of the 

physical realizability condition on K (2.4.11) and the lack of poles 

In the lower half plane (complex) of frequency (Robinson, 1967a). The 

Integration (2.4.13) Is equivalent to evaluating singularities In the 

1/2 

upper half plane only. The Fourier transform of In <|) is given by 


In 4> 


1/2 


a + 2 Z a cos cat 


(2.4.14) 


where the symmetric nature of the spectrum alters the range of summa<* 
tlon. Accordingly, the coefficients, a, are derived from the Inverse 
transform of In 

IT 

“t “ 2 tt I <|j^^^(w) dO) (2.4.15) 

-IT 

From (2.4.15), by forming logarithms of both sides, we have 
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log K(oj) ■ L + Z L cos tot + i £ L sin tot 
° t-1 ^ t-1 ^ 


a + 2 Z cos tot + 10(to) 
° t-1 


(2. A. 16) 


By equating coefficients In cos tot and sin tot we are led to 


L 

o 



(2.4.17) 


L * 2a^ for t > 0 
t t 


(2.4.18) 


The phase 0(to) In (2.4.12) essentially selects a function K with no 
singularities In the lower half plane. The particular phase relation- 
ship is given by 


00 

0(to) « Z L sin tot 
t-1 

00 

« 2 Z sin tot 
t-1 •= 


(2.4.19) 


Equation (2.4.19) is the essential result of this section. The 
analytical problem of spectral factorization Is now complete because 
the phase characteristic that distinguishes the kernel from any other 
function with spectrum (() can be computed In terms of the spectrum Itself. 
This result is clear if we substitute (2.4.15) Into (2.4.19) to produce 


IT 


2 f 1/2 

— Z sin (ot cos (O-t In ((O-) dco- 

t-1 


0((o) 


(2.4.20) 
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The class of spectra to which such an operation will apply depends on 
whether the spectrum obeys the Paley-Wlener condition 


f Is-Miilau, > _ 


(2.4.21) 


This constraint Is discussed In Appendix B as It relates to numerical 
approximations . 

Katzenelson and Gould (1962, 1965) have described a method of ex- 
tending the spectral factorization method to the evaluation of non-linear 
kernels of a functional representation. Their method Involves the 

successive minimization of error between a sampled realization and the 
ttl 

output from an n order representation In order to deteirmlne the 
optimum n+1 kernel. Several hypotheses Implicit In their approach are 
not valid In the problem here. First, the assiamptlon that there Is a 
freedom to generate an output realization of turbulence at will. Is 

not appropriate. Second, the higher order spectra (bispectra ) 

which are needed In Katzenelson and Gould's method are not yet available 
except In very tentative form (Elderkln et al., 1972). Also there Is no 
guarantee that these higher order spectra will obey the factorable 
properties required by the technique. 

We now move to the Implementation of the mathematical development. 
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3.0 LINEAR REPRESENTATION OF UNI-COMPONENT TURBULENCE 


In this chapter, linear models of a single component of turbulence 
that are consistent with surface layer similarity are considered and 
extended to generalized spectral forms using spectral factorization. 
Because the response functions or kernels are a relatively novel con- 
cept In describing the structure of atmospheric turbulence. It Is 
Interesting to study their form under varying meteorological conditions. 
Properties derived from the response function. Including measures of the 
predictability and memory of the model turbulence also are examined. 

The discussion of this chapter will be limited to the vertical velocity 
component, partly for convenience and partly because of Its Importance 
In aeronautical response problems. 

3.1 Model of the Vertical Velocity 

The laqportance of the vertical velocity spectra, both In modeling 
the response of aircraft to turbulence and In studying the vertical 
flux of momentum and heat near the earth's surface. Is Indicated by 
the availability of empirical estimates of Its spectral form. Busch 

2 

and Panofsky (1968) have approximated the w spectra (normalized by u^ ) 
by a form 


f G (f) 
w 


A f 
w 


1+B f 
w 


5/3 


(3.1.1) 


based on data drawn from several sites (f Is defined In (1.4.9)). They 
note that at low wave numbers, f < 1, their empirical form is an 
Improvement over that suggested by Pasquill and Butler (1964), 
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f G (f) 
w 


A f 
(l+B^ 


Confirmation of the Busch'-Fanofsky form Is supplied by Kalmal et al. 

(1972) based on the Kansas data. However, the estimates for (A^,B^), 

for neutral stability, vary between the estimates of Kalmal (1.0, 1.5) 

and Busch and Panofsky (1.5,2. 7). The variation In these coefficients 

Is indicative of the accuracy that can be expected In estimating 

characteristics associated with spectra, such as variance, length 

scale, dissipation, and response functions. All empirical spectra 

-5/3 

behave asymptotically as f , which Is characteristic of the 
Kolmogorov region (1 « f « f j^q^) • 

Another spectral form (Appendix A) 


f G (f) 
w 


A 

w 


2/3 .2/3 


1+(B 

w 



(3.1.3) 


Is convenient mathematically but does not have the usual front slope 
of +1 found by observations. The convenience of (3.1.3) lies in its 
form after a transformation of variables 


f» 



(3.1.4) 


From the Invariance of energy with a change of variables 


f G 


w 


f» 


G* 

w 


we are led to 


(3.1.5) 
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£' G' - =■ (3.1.6) 

1+(B' f) 
w 

which Is a form familiar In filtering problems (Section 2.4) and which 
has a well-defined kernel for a linear response to white noise given by 


K^(x*) « exp - X* > 0 


= 0 


X* < 0 


(3.1.7) 


Because (3.1.4) Is related to a transformation of the space variables 


X* 



(3.1.8) 


the general form of the response function will be 
K^(x) - exp H(x/B^) 


(3.1.9) 


where the function H carries the effect of the linear Fourier trans- 
formation Involved In factoring (3.1.6) to obtain (3.1.9). The 

coefficient V can be shown to be 
'w 


Y - A * 2/3 -2/3 

w w w 


i- 

3tt 


^2 

w 


(3.1.10) 


The response function of any empirical formulation for w spectra, 
to the extent that it approximates (3.1.3), can be expected to vary 
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/ 

according to (3*1.9) • Therefore it Is useful to examine the basic 

structure of response functions in (3.1.9) for H * 1. 

3.1.1 Analytical Characteristics of Response Structure * The 

form of the approximate analytical kernel is shown in Figure 3. 1, 

both as a function of scale distance and of stability, where the 

1/2 

parameter C is defined in (1.4.6). The variation of y (C) and 

w 

B^(C) is based on the empirical relationships of Appendix A and is 
given in Figure 3.2. The general feature of the solution for is a 
mono tonic decrease of the response with distance (in the direction 
from which the turbulence is advected). For a given longitudinal 
separation x, as the height (and scale length, 1) Increase, x • x/i 
decreases. From Figure 3.1, a decrease in x is equivalent to an in- 
creased response. The kernel for separations less than i2'(x < 1) 
decrease faster than exp(- x) but decreases less rapidly than the 

response of a simple linear oscillation system for x > 1. That is, 

~2/3 

the approximate analytical solution, exp(-x ), indicates a decrease 
in the response for small lags (relative to i/u^) , or equivalently, 
indicates that the filter will give less weight to the more Immediate 
past. On the other hand, the response for large lags will be greater 
than that for the common first order linear model. 

The effects of stratification on the kernel, also given in 
Figure 3.1, are two-fold. First, the initial response (x » 0) varies 
with stability, and is a minimum in neutral conditions. Second, the 
rate of decay of response decreases with decreasing stability. The 

Z2 

mlnimtim initial response Is a reflection of the minimum in w at 
C « 0 (Figure A. 3). The rate of decrease of is determined by 
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For Increasing instability (Figure 3.2), Increases as more energy 
Is Introduced at larger f (Lumley and Panofsky, Chapter 5). Therefore, 
from (3.1.9), the response at a given x Increases. 

Another concept which can be demonstrated for the simple 
analytical approximation Is the memory of the system. Intuitively, 
memory may be considered as the Integrated effect of past stimuli. For 
convenience, it Is desirable to compare the memory of the process 
representing turbulence with that of a simple linear first order 
process with the same variance. Memory Is defined tentatively as 


Mem » ■r K^(y) dy} / {K^(0) r exp(- y)dy} 

- 


T^- 1“ 

IT(0) J- 


(3.1.12) 


(y) dy 


It is noted that this definition is only useful if K^(0) ^ 0. For the 
kernel (3.1.9), Mem - 1.33, which Indicates a net Increase In memory 
of about 0.33 relative to a simple first order process. For a 
slightly different version of the memory concept, given by 


Mem(x) 


Im. J 


(y) dy 


K‘(0) 


(3.1.13) 


the memory as a function of distance or time from a stimulus Is Itself 
a function of scale. From Figure 3.1, ^m(x) <1, for x > 1. Therefore 
the Increased total memory which is greater than unity (Mem(<») « 1.33), 
results from the large scale structure of K^. 
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This completes the discussion of the general properties of the 
approximate analytical structure of K^. It remains to contrast this 
Intuitive and preliminary discussion later with more exact numerical 
solutions for a variety of empirical spectra. 

3.1.2 Kernels for Different Empirical Spectra . The spectral 
factorization procedure was applied numerically to some of the 
empirical spectral representations of Busch and Panofsky, Kalmal 
et al., and Pasqulll and Butler as well as the spectrum discussed In 
Appendix A and Section 3.1.1. The objective was to determine the 
range of response estimates which could be expected from variations 
In empirical representations of the w spectra. This variability 
provides a realistic estimate of accuracy against which to contrast 
other sources of variability » such as arise In the parameterization of 
the thermal stability. 

The kernel obtained by factorization of some of the empirical 

spectra are given in Figure 3.3. Also shown are the results for 

arising from the common spectral form used by Busch and Panofsky and 

Kalmal et al. , but with A^ and B^ altered for compatablllty with the 

variance and Kolmogorov range structure (Appendix A) . The final kernel 

plotted In Figure 3.3, and termed "model", corresponds to the spectral 
-1/3 4/3 -1 

form f ^ * with coefficients chosen for compatablllty 

with the variance and Inertial range structure. 

The response function, as expected. Is monotonlc and similar to 
~ 2/3 

the basic exp(-* x ) form of Section 3.1.1. It Is noted that estimates 
of K^ using the Busch-Panofsky formulation differ significantly with 
Increasing scale from either of the formulations based on the Kansas 
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data. The response based on the Kansas data for x > 1 Is not as large 
as the kernel derived from the Busch^-Panofsky data. The variance 
compatible spectrum, based on the common mathematical form used by 
both Kalmal et al. and Busch and Panofsky, results In a response 
structure markedly different from the Busch-Panofsky form alone. It 
Is concluded that the normalized variance characteristics of the data 
set drawn from the Kansas experiment and that used by Busch and 
Panofsky differ significantly. 

The underlying reasons for this disparity are not clear, but may 
be attributed to some degree to the larger roughness characteristic 
of the Busch-Panofsky data set, or perhaps a difference In similarity 
Involving the average structure of the large scale flows (Kerman, 
1974b). Whatever the cause of the disparity In the form of K^, the 
results Indicate that a significant difference exists In the response 
representation at scales, x > 1, resulting from various experiments. 
From Figure 3.1, the estimated errors between the functional forms 
are about equivalent to an error of + 0.25 In an estimate of The 
numerical estimate of the kernel corresponding to the model spectrum 
underestimates the small scale response and overestimates the large 
scale response. An examination of the different spectra factored to 
produce the response estimates (Figure 3.4) reveals the close rela- 
tionship that exists between the relative distributions of variance 
of the spectra and the relative response structures. For spectra 
with additional variance at scales, f < 1, (for example, the analytical 
model spectrum) the result Is an Increased response at scales, x > 1, 


and vice versa. 



43 


The objective of this section was to compare and contrast filters 
arising from various empirical spectra. In summary. It Is concluded 
that perceptible differences In the response structure occur according 
to the empirical representation of the spectra. These differences In 
turn are related to the relative spectral content between empirical 
formulation both for the large and small scale regions. 

The response structure Is next studied as a function of stratlfl'-* 
cation. 

3.1.3 Stratification Effects on Kernels . As discussed in 
Section 3.1.1, the function, K^, could be described by two of Its 
characteristics — Its Initial response at x * 0, and Its integral, or 
memory. The function K^(0) and Its memory Is displayed in Figure 3.5 
as a function of for the particular model described in Appendix A. 
The calculations In the spectral factorization are performed with 128 
points In the Fast Fourier Transform (FFT) algorithm. The response 
at X * 0 differs from that of Figure 3.2 because of the approximations 
made near f (Appendix B) . 

In Figure 3.5, the Initial response Increases with |c| and the 
memory increases monotonlcally with Increasing Instability. The 
response function Implied In Figure 3.5 Is equivalent to that shown 
In Figure 3.1. The response for a given x > 0 Is less for stable 
stratification than for unstable stratification and a minimum for 
neutral stability. The structure of the response function, normalized 
by Its Initial value, also follows from qualitative consideration of 
the change In the spectra with stratification. While the variance, 

Z2 

w (which In (3.1.12) determines the Initial response) Increases In 
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f 


Figure 3.5. Initial response and memory for model spectrum. 
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both stratifications, the spectral bandwidth, B^, (which determines 

the decay rate of response) Increases monotonlcally with decreasing 

stability. Accordingly, the excitement of more large-scale energy 

results In an Increased response at scales x > 1. 

It follows from the scaling of the spectra that the response 

for ah arbitrary stability, scaled by Its Initial response. Is only a 

function of x/B^. For a constant flux layer in which the Mbnln- 

Oboukhov length, L, Is also constant with height, one may equate 

changes In CC** z/L) with changes In height. Therefore, for a given 

X, x/B^ will decrease with height, both because x(* x/il) decreases 

with height and because B ^ decreases with height (Figure 3.2). 

w 

Therefore, the response will Increase as a result of an Increase In 

ZT -- 

w with height and a decrease In x/B^ with height. Under unstable 

conditions, Increases approximately linearly with height, as does 
1/2 

,80 that from (3.1.9), for a given x, 

K^(z) -- z exp(- z (3.1.16) 

1/2 

Under stable conditions, B^ Is approximately constant, while 
again varies linearly with height, so that for a given x, 

K^(z) -w z exp(- z (3.1.17) 

1/2 

For neutral conditions, both y^ and B^ are constants In height, 
and the response function has a form 


K^(z) - exp(- z”^^^) 


(3.1.18) 
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For large enough z, aseumlng Kz/H is still nearly unity, the response, 
for a given x, will vary linearly with height In stratified cases and 
approach a constant In neutral cases. 

It Is concluded that the response structure varies appreciably In 
various stratifications. The results for the Initial response and 
memory, as a function of are Intuitively consistent with qualitative 
discussions of the spectrum. It Is now useful to consider another 
effect on the simulated turbulent process — that of viscosity for low 
Reynolds number. 

3.1.4 Viscosity Effects on Kernels. The empirical spectra dis- 
cussed in the previous sections, were obtained In flows whose Reynolds 
number was sufficiently large that the viscous range was well removed 
from the energy-containing eddies. Accordingly, the models discussed 
there should probably be called 'Inertial* but we will refer to them as 
Invlscld for mnemonic and comparative purposes. It Is useful to 
consider modifications of the Reynolds number (Re) criterion In order 
to study the response structure In the presence of viscosity. This 
problem is not germane to the usual application of filters which 
simulate the energy containing scales. Rather It Is preparatory for 
later attempts (Chapter 5.0) to simulate the derivative structure of 

the small scale region (f «« f „ ) . 

KOL 

The response functions may be considered as the velocity field 
that would be produced by a single liiq>ul8e (Section 2.4)). For 
Invlscld flows, the response to the stimulus Is Immediate, giving a 
discontinuity at the time of the Impulse. In a viscous flow, the 
formation of such Infinite curvatures Is liiq>osslble, requiring that 
response functions rise smoothly to a maximum value. 


I 
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The viscous adjustment to the three dimensional Kolmogorov 

i ■ r 

Spectrum of the Inertial subrange Is Illustrated by the model spectrum 
(Pao, 1965) 

E(k) « exp [- y a(kn"^)^^^] (3.1.19) 

where r| ^ Is the Kolmogorov wave number, (e V ^)^^^. The corresponding 
one~dlmenslonal spectrum for the Isotropic region (Batchelor, 1953, 
p. 50) Is 

4>33(ki) - J^(kj^) (3.1.20) 

where 


J (kj 
w 1' 


II exp[- |a(kj^/C)*^^]de 


(3.1.21) 


and 


k^ - kj^n 


(3.1.22) 


In (3.1.20) 


a 


24 
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a 


(3.1.23) 


where a Is the three-dimensional Kolmogorov constant (-1.5)). A model 
of the spectrum over the entire scale range from energy-containing to 
dissipative Is therefore given by 
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f G (f) 
w 


A (j) f 
w e 


1 + (B f) 
w 


2/3 

4/3 '^w 




(3.1.24) 


where (Appendix A) 


KOL 


1/4 


(3.1.25) 


For a given C» and fixed <j) , the effective cut-off wave number, f » 

£ KOL 

3/4 

of (3.1.24) varies as Re 

The method of digital spectral factorization was applied to 

(3.1.24) for C ^ 0, and several ranges of Re. The response for a given 

value of Re rises from an Initial zero value, overshoots the Invlscld 

case, reaches a maximum, , at a distance x from the origin, 

’ * max’ max ® * 

and then settles down to the Invlscld solution. The variation of 

max 

and x^^^ are given In Figure 3.6 as a function of the Reynolds number. 
Interpreting the figure, we can see that for a decrease In viscosity, 
the maximum response Increases and the displacement of the location of 

3 

the maximum response from the origin Is reduced. For Re * 10 , the 
resulting form of closely resembles the Invlscld result, and the 
maximum response of this viscid case and Initial response of the 
Invlscld case are approximately equal. The location of the maximum 

3 

response for Re * 10 occurs within a distance, 1/100, from the origin, 
or Impulse point. 

2 

However, for low Reynolds numbers, say. Re £ 10 , the effective 

loss of variance due to the viscous spectral cut-off (by J^) Is not 

reflected In the coefficients A^ and B^. Accordingly the computations 

displayed In Figure 3.6 for R^ are underestimates for the lower 

max 

Reynolds number range. 



MAX 
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Another feature of the viscous response structure (not shown here) 
Is the convergence of each kernel. Irrespective of Re, to the 
corresponding Invlscld kernel, for sufficiently large x. The con- 
vergence occurs at progressively smaller x with Increasing values of Re. 
This result Is In keeping with the Intuitive notion that a decreased 
viscosity Is felt at progressively smaller scales where large gradients 
are possible. Equivalently, as the range of unaffected scales extends 
to larger f, the response over a wider scale range becomes Indistin- 
guishable from the Invlscld result. 

The concept of memory (3.1.14) does not apply for a viscous model 

because K^(0) - 0. Therefore, It Is necessary to consider another 

3I/2 

ad hoc normalization, w , rather than the equivalent response of the 
common first order linear model. The memory Is redefined as 


Mem 


0 


K (x) dx / { 


J 

0 


[K^(x)]^ dx}^'^^ 


(3.1.26) 


Estimates of memory by (3.1.26) are given also In Figure 3.6. Varla- 

1 3 

tlon of Re from 10 to 10 results In a 20 percent decrease In the 
memory. This result Is compatible with the Intuitive concept of 
decreased memory with Increased turbulent scrambling as the viscosity 
Is reduced. Another aspect of the memory structure of the simulated 
process Is Its predictability based on Its past history. This aspect 
Is examined next. 

3.1.5 Predictive Structure of the Model . In the development of 
control systems It Is advantageous to be able to predict the turbulent 
velocity field at some future time, on the basis of past observations. 
In order to be able to apply some properties of linear stochastic 
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processes to the linear, Gaussian model of the vertical velocity 
component. It Is necessary to examine another property of response 
functions (Section 2.1). An Inverse linear functional, K Is 
defined by 


j K(tj^) k“^ (T-tj^) dtj^ - 6(t) (3.1.27) 

A Fourier transformation of (3.35) produces an equivalent definition 

/V /S - 

K(w) k“‘‘-(w) - 1 (3.1.28) 


Accordingly, the Inverse linear functional Is derivable from the kernel, 
K, by the method Implied In (3.1.27) or (3.1.28). For a perfect system, 
without noise, the Inverse functional generates a white, Gaussian 
process, from a shaped spectral process, w. In the manner 


C(t) 


j k“^(T) w(t-T) dT 


(3.1.29) 


Let us consider a prediction of w In terms of Its filtered past. The 
filter Is determined so as to minimize the least squares error between 
the prediction and verification (Robinson, 1967b). The linearly 
predicted value, w^(t*Hx), at a time a In the future Is given by 

Wp(t-Kx) ■ I M(x;a) w(t-x) dx (3.1.30) 

For a linear W-H representation of w (Robinson, 1967b), the prediction 
kernel, M, Is given In terms of K and K ^ by the expression 
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M(T;a) - j K(p+a) k“^(t-p) dp (3.1.31) 

From (3.1.27) or (3.1.28), equation (3.1.31) can be considered as a 
relationship between M and K. 

The inverse filter K ^ computed from the kernel of the spectral 
model of Appendix A, under neutral, Invlscid conditions, is given in 
Figure 3.7. The physical effect of K ^ is to filter a correlated time 
series to produce a white noise process. In Figure 3.7, this decoupling 
of the time series is accomplished by the alternate oscillating weights 
near x « 0 in what amounts to a shredding action. A measure of the 
effect of an inverse filter therefore lies in the difference K ^(0) - 
K ^(Ax), where Ax is the resolution for the white noise process which 
will be generated. The larger the difference, the more the necessary 
shredding action to destroy the turbulent correlations. Accordingly, 
for situations with different spectral bandwidth in different thermal 
stabilities the oscillations in K ^ near x = 0 will Increased for de- 
creased stability. 

The results for the prediction kernel, M, are given in Figure 3.8 
for several values of aAx. The most distinctive feature of the 
structure of M is the very rapid decrease in predictive weighting for 
even x < 0.25. The implication is that the best estimate (in the least 
squares sense) of w at a distance oAx ahead is given effectively as a 
multiple of its present value. The error of a prediction aAx units 
ahead is given by 
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S^(aAx) ■ |w(x + aAx) - w^Cx + oAx|^ 


-J 


1 - I M(rAx, aAx) w^((a+r) Ax) d(rAx) 


(3.1.32) 


2 

Intuitively, the rate of growth of error E , with distance Into the 

future. Is a measure of the predictability of the turbulent process 

simulated by the linear representation. 

The error of prediction of the model of Appendix A Is given In 

Figure 3.9. The deterioration of the prediction at even short distances 

~ 2 

Is apparent. For example, at aAx = 0.5, E (0.5) - 0.6, or E - 0.8. 

That Is, at a distance of about f>/2, the root mean square error In 

-jl/2 

estimating the vertical velocity will be about 80 percent of w 

For comparison, a test was conducted of the common first order linear 

process, with comparable bandwidth. The results are shown also In 

2 

Figure 3.9. A comparison of the mean square error E of the turbulent 

-5/3 -2 

spectral model (f ) and the common first order model (f ) Indicates 

a modest Improvement In predictability at scales comparable to It 

Is concluded on the basis of this comparison that modeled turbulence 

Is the more predictable process. This result Is also In agreement with 

the discussions In Section 3.1.1 of the memory of the simulated 

turbulence. 

The basis of the previous discussions Is the response function 
because from It can be derived the memory and predictability char- 
acteristics. An empirical formulation for Is next summarized for 
convenience . 

3.1.6 Empirical Formulae for Simulation Model . In Section 3.1.1, 
the filter for simulation purposes was represented In the form 
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K^(x) - exp-(x/B^)^/^ H(x) (3.1.3A) 

For the purpose of application of the model. It Is convenient to sumr’ 

marlze the forms Y » B and H. 

'w w 

The factor H(x) was computed from numerical results for for 

various stabilities, C> using the following enq>irlcal formulae for 

and B 
w 


- 0.75 (1 + 0.75lc|) (3.1.35) 

* 0.7 (1 + 0.75 C + 3.0 C^) (3.1.36) 


over the range - 1.5 < C < 0.5. The results from H(x) for the 
extremes of the stability range, are shown in Figure 3.10. Apparently, 
H is only a weak function of the bandwidth, B^, of the spectrum used. 
Accordingly, stability effects are ignored in approximately H. The 
form chosen to represent H empirically is given by 


H(x) - 0.5 (1 + exp[- 2.5 x^^^]) 


(3.1.37) 


The specification (3.1.34) of the first order kernel of the linear, 
Gaussian vertical velocity model is now complete. The application of 
these formulae requires establishing estimates of u^ and T^ for a given 
height, roughness and geostrophic wind speed to denormallze the tabu- 
lated functions. An exan^le using the resistance law formulations 
is given by Kerman (1974a) . 
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3.2 Model of the Vertical Velocity Detlvatlve 

In some applications, such as the simulation of the diffusion of 
passive airborne contaminants or the testing of Instrument systems to 
measure turbulence. It may be advantageous to have a simulation of the 
derivative of a velocity conq;>onent. However, the derivative process 
has several distinctive properties which make simulation more difficult. 
First, the maximum spectral content of the derivative of a turbulent 
velocity component lies In the fine scale, viscous sub-range. Second, 
from observations, (for example, Tennekes and Hyngaard, 1972), It 
Is known that the actual turbulent process of the derivative Is dis- 
tinctly non-Gausslan. In this Section the properties of a Gaussian 
derivative process are studied. A somewhat more realistic non-Gausslan 
model Is presented In Chapter 5.0. 

Let the derivative of a linear process, w^ be represented In terms 
of a white, Gaussian process, by 


3w(t) 

at 


j (t-T) ^(T) dT 
0 


(3.2.1) 


A relationship exists between this representation (3.2.1) and the 
representation of the velocity, given by 


w(t) 


I K^(t-T) C(T) dT 
0 


(3.2.2) 


which, by application of Leibnitz's rule, can be shown to be 


D^(t) - K^(t) 


(3.2.3) 
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for K^(0) * 0. Therefore* It Is equivalent to compute directly by 
the spectral factorization arising from the representation (3.2.1) or 
by differentiation In (3.2.3). From the discussion for the spectrum 
near fhe viscous sub^range (Section 3.1.4) the spectrum for the 
derivative* Is given by 

( 3 . 2 . 4 ) 

There Is no need to Include the parameterization of the low wave number 
range In the spectrum of dw/dt because the spectral contributions to 
the derivative are negllble near f » 1. However, the effect of stability 
Is retained In (|) and the viscous effects are parameterized In > 
for a given stability. 

The response function for dw/3t for neutral stability and Re * 10* 

Is given In Figure 3.11 as well as the corresponding kernel for the w 

1 

process. The maximum response* D ^^* In the derivative process occurs 
at the location of the maximum positive derivative In K^. The cor- 
respondence of and through (3.2.3) Is shown clearly. The 
extensive region of negative response In corresponds to the 
monotonlcally decreasing form of over the same range. The resulting 
estimate of the memory (3.1.26) of the derivative process Is slgnlfl~ 
cantly smaller than the memory of the velocity process. This Is 
supported by observations that realizations of the derivative of a 
stochastic process, are more disorganized than the velocity field. 

In Figure 3.12* the maximum response of the derivative process In- 
creases with Increasing Reynolds number. The location* x ^^^ * 


of the 
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maximum response converges to the origin in the same manner as 

max 

approaches the Invlscld limit for large Re In the velocity process. 
However the negative response also Increases In width and magnitude 
with Increasing Re so that the memory continues to decrease. Therefore, 
different time or space series of a realization of dw/dt would display 
Increasing disorder with Increasing Re. Qualitatively, however, such 
realizations would not display the patchy nature often observed In real 
turbulent realizations (Stewart, 1969). An actual turbulent derivative 
signal would Include more large negative excursions (rapid decelera~ 
tlons) than large positive excursions and have more large excursions of 
any sign than a purely linear Gaussian process. An analysis of these 
features is given In Chapter 5.0. 



4.0 LINEAR REPRESENTATION OF MULTI-COMPONENT TURBULENCE 


The development of the previous chapter assumed that a velocity 
or derivative component was uncorrelated with any other component. But 
requirements exist, or may be foreseen for simulating the longitudinal 
velocity, u, the vertical velocity, w, and the buoyancy, 0, components 
of a turbulent flow near the ground. For example, the requirement for 
a multi-component model of turbulence for VSTOL response problems has 
been outlined by several authors (Case, 1968; Skelton, 1968; Houboult, 
1973). 

The method which Is used to calculate multi-component filters 
Involves the reduction of the spectral matrix to an equivalent series 
of single component spectra, with subsequent spectral factorization. 

The reduction process Is discussed next. 

4 . 1 Model Development 

Consider a multi-component, linear representation In terms of a 
white, Gaussian process for the longitudinal and vertical velocity 
components and buoyancy. In the direction of the mean wind speed 

u^(x) “ I K^j(x*) C^(x-x') dx’ (4.1.1) 

For later consistency, It Is convenient to define a velocity, u^, from 
the buoyancy fluctuation 
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80 that Is understood to be the array (u, w, 

represents the response of the 1^^ output velocity component to the 
white, Gaussian Input conq}onent. It Is not necessary to Include 
the lateral velocity component, v, because It Is uncorrelated with all 
other variables and can be simulated Independently. 

A Fourier transformation of (4.1.1), results In the expression 

/S 

5^(f) - Ky(f) Cj(f) (4.1.3) 

and the spectral matrix la given by 

^mn m n 

TV /V 

- K - 

mj Jn 

✓V 

Because the process, (!) , is Hermit Ian (K. . 

mn ij 

equivalent diagonal process such that 

4» - n (4.1.5) 

where S) Is diagonal. 4 Is the eigenvector metrlx of (|>. (The symbol t 
represents the adjoint, or transposed complex conjugate of a matrix.) 
Because the dlagonallzatlon procedure Is a linear operation, equation 
(4.1.5) applies at each scale Independently. The diagonal elements of 
Q((«)) represent an uncorrelated process equivalent to a single-component 
spectrum. Therefore, we consider next the.,^spectral factorization of 


(4.1.4) 


), there exists an 




The kernel K 


ij 
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each elgenspectrum (or diagonal element of f2). The result of the 
operation on uncorrelated diagonal elements Is Itself a diagonal matrix 
(of spectral factors), say, X, where 

n - X x'*’ (4.1.6) 


From (4.1.4) and (4.1.5) 

= (♦ K) (<|) K) 


(4.1.7) 


so that 


Xy - K 


(4.1.8) 


where y Is an arbitrary unitary matrix. After essentially a trial and 
error analysis, It was found that the equivalence 

y - * (4.1.9) 

preserved the minimum phase characteristics of K. Consequently, from 
(4.1.8) and (4.1.9) 

X - * K (4.1.10) 

or, the eigenvectors, diagonalize both And K. 


This conq>letes the mathematical development of the multl'-component 
model. In summary, the procedure Is first to reduce a matrix of spectra 




Figure 4.1. Spectra, co-spectrum and eigen-spectra of two 
component model. 
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on the other hand, underestimates the w spectrum In the large scale 
region but converges to the u spectrum In the Inertial sub-range. 

The response functions corresponding to the elgenspectrum (termed 
elgenresponse for convenience) are similar In form to the single 
component model. The multi-component functions derived from the elgen- 
response functions by rotation using the eigenvectors (4.1.10) are 
displayed In Figure 4.2 The response, ^ first Gaussian 

Input Is similar to the first order response function, . The kernel 

u 

has a larger response at small scales (x > 0.5) than a single 
component model, but a smaller response at larger scales (ic > 0.5). 

The effect of a downward momentum flux on the vertical component Is 
to decrease the response for all scales. The response ^^31^ 

Is negligible for small scales and becomes approximately constant 

~ 11 
for X > 0,5 Although both the self-responses, K and ^ 331 

reduced by the presence of the stress, the total response for 

simultaneous, equal and opposite Impulses In the Input channels will 

be larger than If the u and w consonants were uncorrelated. 

The predictability of the multi-component model was also studied. 

The mean-square error for the single components and the mean square 

error of w for a multi-component simulation are plotted In Figure 4.3 

The errors of prediction are significantly different between the u and 

w components In single component models. The w component Is Inherently 

less predictable and the lnsosltlon of a cross-correlation makes only 

a minor change to 2^ , Equivalently, the prediction functions 

w Xj 

are only slightly different than the single component estimates. 

It Is concluded that the response for x 0.1 for a multi-component 
mode in neutral stratification may be increased by the presence of a 
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downward momentum flux. However) predictions based on u and w slmul** 
taneously show only marginal Improvement over single component models. 

4.2.2 Model In Non~Neutral Stratification . The basic u--w model of 
Section 4.2.1 was enlarged to Include u-0 and w-0 correlations. Co- 
spectra Involving the lateral, v, component were neglected but the v 
spectrum was Included. The structure of characteristics of the response, 
memory, and predictability for a multi-component model under various 
stratifications, are discussed next. 

The distribution of Initial response (x = 0) of the eigen-response 
functions for various stabilities (Figure 4.4) is similar to the single 
component response functions. The, ranking and indexing of Is by 
magnitude. The minimum response of the eigenstructure Is the same as 
for the single component w model In neutral stratification. Whereas 
the initial responses of the first and second orders are identical, 
their memories (Figure 4.5) differ. Accordingly, the order 1 and 2 
eigenfunctions represent two distinct processes but with the same 
Initial response. The elgenmemorles Increase generally with increasing 
instability — a property shared by single component models. 

The transformation from the eigen-response to the multi- component 
response structure by the use of the eigenvectors results In the same 
qualitative picture as outlined for the eigen- and single-component 
structure. As before, initial responses are greater in non-neutral 
stratification with all variables, and the memories increasing with 
increasing Instability. The Initial response structure of the v 
component is identical to that of the w component because of their 
Identical inertial sub-range structure. 
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Figure 4.4. Initial eigen-response for four component model as 
a function of thermal stability. 
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Figure 4.5, Eigen-memories for four component model as a 
function of thermal stability. 
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Figure 4.6. Extrema o£ response of stress and heat fluxes as a 
function of thermal stability. 
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Figure 4.7. Mean square error of prediction of vertical 
velocity for various stabilities. 
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The cross-response K (i j ) Is a measure of the coupling between 
the components, u^ and u ^ . The major cross-response functions each 
have zero Initial response, similar to the result of the previous 
sectloit. Another feature of the absolute value of the cross-response 
functions Is a maximum value obtained at a scale comparable to x * 1 . 

These extreme values for ^^ 3 * ^^ 3^9 ^^30 plotted in Figure 4.6 

as a function of stability. The response function, represents the 
response of u^^^ to ^nd therefore the coupling be- 

tween u and w. From Figure .4.6, the coupling between u and w decreases 
with Increasing Instability. This result Is In accordance with the 
results of Wyngaard et al, (1971b) who deduced that the approach to 
free convection Implies a loss of preferred longitudinal, or x, direction. 
Their argument was based on the negligible value of u 6 for ^ > 1. On 
the other hand, the absolute maximum values of the response functions, 
and seen In Figure 4.6, Increase away from neutral 

stability. This represents an increase in coupling in stratified flows, 
initially between w and the buoyancy, and subsequently between u and 6 . 

The effect of thermal stratification on the predictability of the 
multi-component flow was also studied. The mean square error of the w 
component as a function of the prediction distance aAx and several values 
of stability Is shown In Figure 4.7. The predictability of w In neutral 
cases was discussed In the previous section. The predictability of w 
Increases with thermal Instability and decreases with stable conditions. 

It Is concluded that the predictability Is Increased or decreased by 
the Increase or decrease of low frequency content of the flow under 


various stabilities. 



This con^letes our discussion of linear » Gaussian models and we 
move to a discussion of several non-llnear and non-Gausslan models. 
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5.0 NON-LINEAR REPRESENTATION OF TURBULENCE 

Observations o£ turbulence In the atmospheric boundary layer, as 
reviewed In Section 2.2 and Appendix A, have established that the 
velocity structure of the energy containing region of the spectrum Is 
nearly Gaussian. This fact Is exploited In this chapter to build a 
weakly non-llnear model to simulate the low order moments and spectrum 
of surface layer turbulence. The development In terms of W-H functionals 
allows for systematically Incorporating some characteristics of tur- 
bulence, such as Its so-called patchiness. Previous simulations of 
the non-Gausslan structure have relied on ^ hoc methods (Dutton, 1970; 
Reeves, 1969). 

The extension of filtering methods to non-linear models is not 
without Its difficulties, particularly with respect to the lEq>llcatlon8 
about the dynamics of turbulence. Models of turbulence In the energy 
containing region and of the derivative structure In the fine scale 
spectral region are developed and studied In this chapter. The 
derivative model Is shown to simulate the transfer of energy towards the 

viscous subrange In^a manner similar to a theoretical model of Pao 

(1965, 1968). 

5.1 Velocity Model 

In Section 2.3, the equivalence between the statistical structure 
of turbulence In the form of moments of the probability density 

function, and In the form of kernels of a functional representation 

was examined briefly. Let us pursue this equivalence somewhat further. 
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The W-*H functional expansion of a (one-dlmenslonal) representation of 
velocity Is given by 

u(x) « Hj^Cx - x^) dx^ + *2^ 

H2(x - x^, X - x^) dx^ dx2 + *2* *3^ (5.1.1) 

H^Cx - x^, X - x^, X - x^) dx^ dx2 dx^ + ... 

A hierarchy of moments follows from (5.1.1) by forming successive 
products between u(x) and u(x + Ax) and averaging In x. The result Is 
a set of simultaneous Integral equations In the kernels, which, for 
Ax « 0, Is given symbolically by 


u - 0 


u^ - j • j + 2 *1 '''*•* (5.1.2) 

u^ - 3 I • I 


M 


2 3 3 

+ . . . 


(5.1.3) 


u^ - J • I + 4 I • I ^ j ’ 1 


I • j kVk^k^ '*‘1*1 


(5.1.4) 


L 2 2 3 

K K + . . . . 
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In the derivation of (5.1.2) - (5.1.4), the orthogaonallty of the 
Hermlte polynomials, H^, with respect to the p.d.f., (Barrett, 1963) 
Implied In the averaging, has reduced the many cro8s-*product8 among 
the polynomials to simple Integer coefficients which Indicate the 
miiltlpllclty of products of the Dirac 6~functlon. A corresponding 
reduction In the order of the Integrations Is also Implied In the 
symbolism 

The hierarchy of moments Is Infinite both In the number of moments 
and the number of terms In each series. However, as observations only 

~4 

exist up to u (Appendix A) for^the atmospheric boundary layer, the 
problem reduces to a finite sub-p.roblem. Consequences of this method 
of closure will be examined In analyzing the results of the model. 

The formulation of the resulting finite set of coupled Integral 

3 

equations (truncated at K ) Is given by (5.1.2) - (5.1.4). These 
equations are not easily amenable to solution without further slmpll- 
cation. A convenient method of approximation which overcomes many of 
the computational difficulties associated with the Integral equation 
structure Is an assumption that higher order kernels are multiples of 
the linear kernel, say 


I K^K^. 


(x^. . . .Xj) 


TT 

k=l 




(5.1.5) 


We shall refer to this approximation as the method of separable kernels. 
Equation (5.1.5) Is analogous to the turbulent closure schemes In 
which higher order statistics are expressed In terms of lower order 
statistics. For example. In the classical closure scheme of most 



82 


boundary layer models, the stress Is expressed In terms of the shear 
of the mean flow. Similar arguments are proposed In the quasi-normal 
models of Donaldson (1972) and Deardorff (1972b) for the vertical energy 
flux (a triple velocity product) In terms of the local energy (variance) 
structure. 

Physically, the separability hypothesis (5.1.5), when applied, say, 

2 

for K , implies that the conditional response at the present Instant 
for an Impulse at x^, given a previous Impulse at X 2 » Is given by the 
response at x^ weighted by a multiple, A^, of the response at X 2 * 
However, the conditioning of the response at x^ by a multiple of the 
response at a time, or separation ^ 2 “^^ distant, is physically un- 
realistic. Instead, the weighting of the Impulse as x^ would depend 
more likely on the Integrated history of the response from X£ up to 

2 

x^. Accordingly, a more physically consistent closure scheme for K 
would be 




fX 


1 „1 


K (x) dx 


(5.1.6) 


This representation will be discussed again In terms of a separable 
kernel model for the derivative. There, the concepts leading to 
(5.1.6) are shown to be somewhat easier to specify In a phase space 
representation. Just as the simple Newtonian stress-shear relation- 
ship Is often questioned (Lumley, 1970) on the basis of local 
representations, the overwhelming practicality and reasonable experi- 
mental agreement demand It be retained in lieu of a workable alternative. 
As demonstrated next, the simplicity of local separability (5.1.5) 
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reduces an otherwise unmanageable problem to workable proportions. At 
the same time. It Is possible by this method to achieve a reasonable 
simulation of the moment and spectral structure of the turbulence of 
the atmospheric boundary layer. !The separability conditions 


- A2 kV (5.1.7) 

K^kV" (5.1.8)* 


are substituted In the truncated moment expansions (5.1.2) to (5.1.4). 
The result Is 

u^ - + 2A2 + 6Aj (5.1.9) 

u^ - 6A2L^ + 8A2 L® + 36 A2A3L^ + 108 A2A^L® (5.1.10) 

u^ - 3L^ + 24 A^L^ + 60 A^L^ + 60 A^L® (5.1.11) 


where 


- I [K^(x)]^ dx, - (L^)^, (5.1.12) 

The utility of the assumption of local separability Is now clear. The 

problem has been reduced from a problem In simultaneous Integral 

equations to a problem In simultaneous algebraic equations. Further, 

2 

the e:q>anslon Is In terms of L , which from (5.1.9) can be seen to be 
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the linear contribution to the variance. Therefore, for a weakly non-- 
linear functional representatln of a weakly non-Gausslan process, the 
moment expansions are perturbations about the linear, Gaussian process. 
Further, It Is reasonable to describe such a representation as a 
quasl-llnear model. 

A similar expansion for the spectrum arises from the separability 
conditions. 


lK^(f - dfj^ 


+ II |K^(f2> 


|K^(f - - f^)!^ dfj^ dfj 


(5.1.13) 


The spectrum is then also expanded In the linear contribution to the 
spectrum, |k^|^. 

The distribution of the moments (variance, skewness and kurtosls) 

of velocity and buoyancy, as a function of thermal stability, are 

described In Appendix A. The data are drawn from the Kansas experiment. 

2 

The solution for L , A^ and A^ from (5.1.9) to (5.1.11) will therefore 
also be a function of stability. For convenience, the model is re- 
stricted to a non-Gausslan extension for the vertical velocity component. 

Equations (5.1.9) to (5.1.11) were solved by Iteration. The results 
for the contributions to the variance, skewness and kurtosls for 

various thermal stabilities are given In Figures 5.1 to 5.3. 

2 

In Figure 5.1, the linear contribution, L , Is the principal source 
of variance. The non-linear contributions (2A^ and 6A^ L^) Increase 



Figure 5.1. Partition of variance of vertical velocity for a 
cubic, separable model. 





Figure 5.2. Partition of skewness of vertical velocity for a 
cubic, separable model. 





87 
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Figure 5.3. Partition o£ kurtosis of vertical velocity for a 
cubic, separable model. 
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Figure 5 
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311 of Spectrum of vertical velocity for a 
separable model. 
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monotonlcally with decreasing thermal stability. The cubic contribution 

to the variance exceeds the quadratic portion. In Figure 5.2, for the 

partitioning of the skewness, the major contributor at all C Is from 
4 6 

the term • The term becomes only an appreciable fraction 

of the total skewness for ^ - 1. However, In the contributions to the 

kurtosls (Figure 5.3), the cubic term, 24A^L^, dominates the quadratic 
2 6 

term 6 OA 2 L . 

In summary, the major contribution to the (non-zero) skewness 
Involves a term linear in » as shown In Figure 5.2, and the deviation 
of the kurtosls from a Gaussian representation is accounted for by a 
term quadratic In A^, as shown In Figure 5.3. However, the major non- 
linear contribution to the variance (Figure 5.1) occurs from the third 
order effect in A^ rather than a quadratic effect in A 2 * Alternatively, 
the distribution of variance Is not monotonic with the order of the 
kernel. This effect has been reported and examined by Crow and Canavan 
(1970) . 

The spectral partition of variance (5.1.13) resulting from an 

2 

Iterative solution using the values of L , A 2 and A^ for neutral 
stability is presented In Figure 5.4. Both the quadratic and cubic 
contributions Increase with decreasing scale. The Increased effect of 
non-llnearltles with Increasing f is qualitatively consistent with the 
results of Stewart (1951) and Frenklel and Klebanoff (1967). From 
dynamical considerations. It Is expected that non-linear effects are 
dominant In the Inertial sub- range, say f 10. However, the largest 
relative contribution to the variance (non-llnear/llnear) at f 10 
Is only about 0.25. Although this ratio Increases with stable 
stratification. Its maximum value under any stratification is only 0.5. 
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Therefore, there are Indications that the dynamical structure of the 
model Is not represented properly. This point Is examined further In 
a later discussion. 

In Chapter 3.0, the concept of memory, for a linear system, was 
defined by 


Mem 


j K^(x)dx/j 
o o 


(5.1.14) 


An alternate concept for memory, based on (5.1.14) Is the limiting 
response as x of the system to a step function at x ■ 0, that Is, 

4 


^(x) - 1 (x > 0) 


(5.1.15) 


- 0 (x < 0) 


(5.1.16) 


11 2 

The normalization by [K (x) ] dx simply redefines the system as one 
having an output with unit variance. It Is natural to extend this 
alternate concept of memory to a non-linear system. The resulting 
definition of memory Is given by 


Mem 


{Mem^ + A2(Mem^)‘ 


+ A3(Heiiij^)^}/u^ 


(5.1.17) 


where 



K^(x)dx 


(5.1.18) 
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In (5.1.17) 


2 2 222 223 

u - + 6A3(L^)^ (5.1.19) 

Equation (5.1.17) is consistent with previous decompositions because 
the total memory Is Itself expanded In terms of Its corresponding linear 
contribution, Mem^. 

The decomposition of the memory for the non-linear representation 

2 

of the w process, using the solution for L , A^ and A^, Is presented In 
Figure 5.5. The linear memory contribution, Mem^ (defined in (5.1.18)) 
Is less than that for a totally linear representation. The loss of 
linear memory follows from Figure 5.4 for the spectral distribution of 
linear variance. Because the ratio of non-linear variance Increases 
with f, the contribution to near ic >■ 0 must be less for a non-llnear 
model than a linear representation. Accordingly, the linear memory Is 
reduced. The linear contribution, Mem^, represents more of the total 
memory In stable stratification than In unstable stratifications. This 
trend Is consistent with previous results for the increasing effect of 
the non-llnear aspects of the functional representation for decreasing 
In Figure 5.5, the relative contributions of variance In the quad- 
ratic and cubic terms Is also reflected in the relative contributions 

of the cubic and quadratic terms of the total memory. 

1 

The response function, K , of a non-llnear representation departs 
moat from the corresponding fimctlon of a linear representation x « 0. 
The non-llnear memory changes are exaggerated by the squaring and cubing 
operations In (5.1.17). For Initial responses greater than unity, (such 
aa occur for C < 0) the result Is an Increase In the contributions to 
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i 

Figure 5.5. Partition of memory of vertical velocity for a 
cubic, separable model and a comparison with 
linear model. 
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the non-Ilnear memory components at the expense of the linear component. 
Conversely, for initial responses less than unity, (as for the stably 
stratified region, C ^ 0), the effect of squaring and cubing the linear 
menoTf contribution la to reduce the non-linear memory contributions 
relative to the linear part. Both of these conclusions are consistent 
with the relative distributions of memory in Figure 5.5. 

An energy cascade by interactions among wave numbers is an inq>ortant 
characteristic of turbulent flows with a large Reynolds number. There- 
fore, it is interesting to estimate the inertial transfer of energy 
arising from the non-linear term of the equation of motion (u3u/dx) for 
the simulated one-dimensional process. First, let us consider the 
Fourier transform of the inertial term. The result is 


Tr(f) 


- fjf {u*(f) 


‘ /\ 

{u(f^) u(f - f^) df^} 


(5.1.20) 


where Im{} represents the imaginary part of a complex argument. For 
the W-H expansion of u, given by (2.3.6), the expected value of 
(5.1.20) becomes 


Tr(f) - - 3 j K^(-p) K^(f,p) dp 

- 2 j ^ K'^(“P) K'*’(p-f) dp 

1 r f -^2 ^2 

+ — J K (p, f-p) J K (-P, r) K (-r, p-f) dr dp 


(5.1.21) 


+ . . . . } 
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Further, the separability condition, (5.1.7) reduces (5.1,21) to 
&(f) - - Im {AjlK^Cf)!^ I |k^(p)|^ dp 

IT ^ 

- ^ I |K^(p)|^ iK^Cf - p)|* dp (5.1.22) 

+ ^ j‘-|K^(p)|^ l^(f - P)|^ I |K^(r)|^ 


dr dp + . . , } 


An Inspection of (5.1.22) reveals that the terms bracketed by {} are 
all real. Therefore, Tr Is given as the Imaginary part of a real num- 
ber, and consequently 


Tr(f) - 0 


(5.1.23) 


Accordingly, there is no transfer of energy at all by the separable 
kernel model of velocity. In Its present form, the quasl-llnear model 
Is dynamically Inconsistent with known characteristics of turbulent 
flow. 

In order to have a non-trlvlal transfer of energy. It Is necessary 
to have an Imaginary part to the terms bracketed by {} In the right 
hand side of (5.1.21). The fact that the separable kernel model has 
only a real part stems from the lack of phasing between the Fourier 
transforms of the kernels. The Interactions between kernels have been 
shown by Crow and Canavan (1970) to be the equivalent mechanism for 
transferlng energy as the Interactions between velocity components at 
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different wave numbers. These authors also demonstrated that the W-H 
functional representation Is handicapped by a need to Include very 
many kernel interactions to simulate the Interactions of u over wave 
number space. Therefore, the quasl-llnear model for velocity shares 
the dynamical shortcomings of W-H functional representations that are 
significantly more sophisticated. 

However, the redeeming feature of the separable kernel method Is 
that It leads to a faithful reproduction of the low order statistical 
structure of observed data. We next examine a model for the derivative 
of the longitudinal velocity component which better represents the 
dynamics of energy transfer. 

5.2 Derivative Model 

In earlier discussions (Section 2.2), it was noted that the small 
scale (f » 1) structure of the derivative of velocity In isotropic 
turbulence was distinctly non-^Gausslan, and that the dynamical process 
at such stales was distinctly non-llnear. The results for the non- 
linear velocity model of the previous section Indicate that a simulation 
of the dynamics of turbulence with a truncated set of kernels, even for 
a weakly non-Gausslan situation, requires a more sophisticated closure 
scheme. We now demonstrate that the method of separable kernels applied 
to a model of the derivative of longitudinal velocity, for a modest 
range of Reynolds number, results In a plausible dynamical analogy. 

Consider the W-H functional expansion of a derivative of the 
longitudinal velocity In the mean wind direction 
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1^- d(x) - I d\ + || H 2 + ... 


(5.2.1) 


The relationship of the moments of (5.2.1) In relationship to the 
kernels, Is similar to that given In (5.1.2) to (5.1.4). In the 
same abbreviated notation, the Interrelationships are 


- I • I + 2 j • I dV + 6 j • I dV + ... (5.2. 

d^ - 3 I • I ^^1*1 


2) 


(5.2.3) 


+ . . . 


d^ - I • j '*’^1*1 '*’^1*1 

dVdV + 1*1 ‘*’1*1 (5.2.4) 


J • I .V 


+ 6 I • ! D*D^D^D^ 


where the orthogonality of the Hermlte polynomials has determined the 
Integer coefficients. Again, In analogy with the non-llnear velocity 
model. It Is convenient to assume that the kernels of the derivative 
process are separable, 

1 ^1 

D*^(x. ,...x ) - A TT D-^(x ) (5.2.5) 

^ ^ J-1 J 

The coefficients correspond to the coefficients of the velocity 
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model, ^le result of substituting (5,2.5) in (5.2.2) to (5.2.4) is a 
set of algebraic equations 

+ 2Aj A^ + 6 A^ + ... (5.2.6) 

d^ - 6 A 2 A^ + BA^ ^ + 36 A^ A^ A^ + 108 A^ A® 

(5.2.7) 

+ . , , 

d* - 3X® + 24 A, X® + 60 A^ X® + 60 A^ X® + . . . (5.2.8) 

2 

where A is the linear contribution to the variance of the derivative 
given by 




(x) ] dx 


(5.2.9) 


The structure of moments for the derivative, up to d , are drawn 
from results reported by Wyngaard and Tennekes (1970). Their results 
for the skewness and kurtosis for the longitudinal velocity component 
form a convenient basis for a truncated representation similar to that 
given in Section 5.1. Based on the log-^normal probability density 
function for a derivative process (Gurvich and Yaglom, 1967) , and 
Kolmogorov's (1962) hypothesis concerning the probabilistic structure 
of the local dissipation, Wyngaard and Tennekes found from some 
empirical data that 


3/2 


3/16 




d®/(d^) 


- - Cl Rj 


(5.2.10) 
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(5.2.11) 

where is Reynolds number based on the Taylor microscale. Examination 
of Figures 5 and 6 of the paper of Wyngaard and Tennekes lead to the 
en^lrlcal estimates 




1/2 


- 0.16 (5.2.12) 

Cj - 0.40 (5.2.13) 


Further, the Taylor microscale, can be related to the generative 
scale length, it, for turbulence In local balance between generation 
and dissipation (Tennekes and Lumley, 1972, p. 67) 




(5.2.14) 


and consequently 


Rj - (15 


(5.2.15) 


In terms of the Reynolds number of the generative region. Re, the 
skewness and kurtosls of the derivative become 


d3/(d2)3/2 


- - 0.21 Re 


3/32 


(5.2.16) 


dV(d^)^ - 0.79 Re^'^^ 


(5.2.17) 
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Wyngaard and Tennekes pointed out that dependence of the skewness 

and kurtosls on the Reynolds number violates the original hypothesis 

of Kolmogorov. However, the dependence on Re Is rather weak because 

of the smallness of the exponents, 3/32 and 1/4. To be consistent with 

the formulation leading to (5.2.16) and (5.2.17), it is necessary to 

re~lnterpret the dissipation in terms of a locally averaged estimate 

(Itolmogorov, 1962; Oboukhov, 1962). 

The estimates (5.2.16) and (5.2.17) for a given value of Re were 

2 

substituted in (5.2.6) to (5.2.8), and the equations solved for X , 
and A^ by iteration. The model must be restricted to neutral stability 
because there is no available information on the distribution of the u- 
derlvative skewness and kurtosls with stability. Table 5.1 summarizes 
the distribution of variance, skewness, and kurtosls for various Re. 

For Re 10 , the iterative method of solution of the non-linear 
algebraic equations, (5.2.6) to (5.2.8), did not converge. The skew- 
ness and kurtosls estimates according to (5.2.16) and (5.2.17) imposed 
by truncated of the moments to a finite number in order to represent a 

3 

process with a large Reynolds number (>10 ) is not uniformly valid. 

Also, Inspection of Table 5.1 indicates that the non-linear contribution 

6 2 

to the kurtosis (24 A^X ) is negative for Re £ 10 . Because only 
positive values of kurtosis, like variance, have physical significance, 
this result is Indicative of an unrealistic interaction among the ker- 

3 

nels remaining after truncation. However, for a range of Re between 10 
4 

and 10 , where the model does not radically depart from a Gaussian model, 

2 

the solution for X and A 2 and A^ appears to be successful. 


An empirical form 



Table 5.1 


Partition of variance, skewness and kurtosis of u - derivative 
Variance Skewness Kurtosis. 


logjO 


2 2 


6A2X* 

SA^X* 

aeAjAgX® 

3X^ 

24A2X* 

eoA^x® 

1 

,958 

.012 

.030 

-.443 

.000 

.191 

2.64 

-1.58 

.342 

2 

.987 

.008 

.004 

-.377 

-.002 

.062 

2.89 

-0.64 

.240 

3 

.966 

.005 

.029 

-.276 

.000 

-.117 

2.71 

1.59 

.131 


IpO 
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A f 
u 


5/3 


- r+iB-f)- 


(5.2.18) 


was chosen for the spectrum of the u *- derivative. L Is the one dimen- 
sional# similarity solution (Fao, 1965) for the viscous subrange, given 
by 


L(f) 


(1 - exp (- I a(f/5)^^^] 


(5.2.19) 


where 


f - f (5.2.20) 

for neutral stability. The coefficients A and B In (5.2.18) were 

u u 

calculated according to the variance and dissipation compatibility 
conditions (Appendix A) . 

The linear contribution to the spectrum of the derivative is given 
In Figure 5.6 for several values of Re. The effect of the viscous 

3 

tall does not become significant In the case of Re *• 10 until f - 350 

which Is considerably smaller than the scales of Interest for the majority 

of engineering applications. However, the maximum values of Re and f 

of the model are underestimates of the respective atmospheric values. In 

Figure 5.6, there Is an Increase In the variance of the derivative at 

smaller scales with an Increase In Re. The Increase In the variance Is 

Itself felt In the Increase In the response at smaller scales (Table 5.2). 

The structure of the response function, D^(x), Is given In Figure 5.7 
2 

for Re • 10 and 5*0. The response functions for a total linear 
representation of du/dx and the linear and non-linear representation of 









I 



X 

Figure 5.7, Comparison of non-linear response of derivative 


of longitudinal velocity component with linear 
response of derivative and with non-linear and 
linear responses of longitudinal velocity 
component . 
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Table 5.2 

Magnitude and location of the maximum response 
and linear contribution to the memory of u - derivative model 


logj^O 


X 

Mem. 

max 

max 

1 

1 

0.25 

0.20 

0.824 

2 

2.8 

0.027 

0.0 


u are also presented. The maximum response for the u - derivative 
process occurs at large separation and Is less than the maxlmtun response 
for the linear process. The reason for changes in the non-linear 
structure relative to the linear were discussed In Section 5.1 In re- 
lation to the non-linear velocity model. The increase of the ratio of 
non-linear to linear variance with decreasing scale, for a given total 
spectral content. Is associated with the reduction In the linear response 
at small separations, say x < 0.5. 

In the previous section, we discussed the lack of dynamical 

consistency for a separable model of the non-linear, non-6ausslan 

representation of velocity. The lack of phase Interactions among the 

^1 

Fourier transform of the response functions, K , resulted In a lack of 
energy transfer Inertlally. Therefore, we examine the Inertial transfer 
properties of a separable W-H representation for the derivative for Its 
phase and energy transfer properties. A aeries of equations for the 
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Interrelationships between the kernels and follow from the W-H 
eiq>anslons for the velocity (5 *1.1) and the derivative (5.2.1) when 
each expansion Is substituted In the Identity 


u(x) 




(5.2.21) 


Terms of the expansion are gathered In like orders of the (orthogonal) 
Hermlte polynomials. For the separable derivative model, the Inter- 
relationships become 


K^(Xj^) - I (y + dy 

1-00 

K^(x .X^) * ^2 1 *2^ 

^ • mOO 


(5.2.22) 


(5.2.23) 


The relationship (5.2.23) between the kernels of the corresponding 
velocity model differs from the model of Section 5.1. In particular 


K^(Xj^,Xj) ^ Aj K^(Xj^) K^Cxj) 


(5.2.24) 


where Is a constant. 

A procedure similar to that used In (5.1.1) and (5.2.1), but using 


Fourier transformed equivalents, results In a sequence of expressions 
''l '^l 

between D and K given by 


1 (fj^) (5.2.25) 

1 (fj^ + fj) K^(f^,f2) - ©^(fj.fj) 


(5.2.26) 
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Substitution of the Fourier transform equivalent of the separability 

^2 

of the derivative process (5.2.5) leads to an expression for D In terms 
‘'I 

of D given by 


D 


2 


(fl.^2^ 


^2 ^^1^ ^^2^ 


(5.2.27) 


Finally, the combination of (5.2.25), (5.2.26) and (5.2.27) results In 
a succinct statement of the kernel Interaction structure of the cor- 
responding velocity model 





^ ^ 1 . h 1 1 

^ (fn) ^ (fo) 
1 . 2 ^ ^ 


(5.5.28) 


From (5.1.21), the energy transfer for the separable derivative model 
becomes 


Tr(f) - ^ Aj Im {1 f|K^(f)|^ [ ^ |K^(p)|^ dp 

(5.2.29) 

^ I p(f - p) |K^(f - p) 1^ |K^(f - p) 1^ dp} 

'^2 

An additional term Involving a triple product of K (and hence ^ 2 ^ 
been neglected In (5.2.29). The expression In (5.2.29) In the brackets 
{} Is purely Imaginary In contrast to (5.1.21). Therefore, the energy 
cascade Is non-trlvlal for a separable derivative model In contrast to 
that for a separable velocity model. For convenience, we simplify 
the notation by Introducing a function for the common term f|K | , 


A(f) - 


(3.2.30) 
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Equation (5.2.29) reduces to 

lF(f) - ^ {A(f) I ' 2f I A(P)A(f-p) dp} (5.2.31) 

Computationally, (5.2.31) Is evaluated from estimates of A based on 
(5.2.30) and (5.2.25). 

Estimates Of f Tr(f) derived from (5.2.31) are displayed In 
2 

Figure 5.8 for Re * 10 . The transfer of energy Is confined to the 
viscous region, with negligible energy transfer In the range 10^<f<10^. 

For comparison, we have Indicated the theoretical prediction of the 
transfer function arising In the model of Pao (1965, 1968) for Isotropic 
turbulence. Pao models the energy transfer as a convergence of spectral 
flux, S(f), given by 

S(f) - exp - 3/2 a (5.2.32) 


where 


Tr(f) * ^ S(f) (5.2.33) 

The curve for Tr shown In Figure 5.8 represents a one~dlmenslonal 
energy transfer estimate derived from (5.2.23) based on the assumption 
of Isotropic energy transfer (Batchelor, 1953, p. 50). The agreement 
Is gratifying although the Pao model contains a wider Inertial sub- 
range of nearly zero energy transfer than does the separable derivative 
model. 



fTr(f) 
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Figure 5.8. Comparison of spectral energy transfer In viscous 
sub-range for derivative model and similarity 
model of Pao. 
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The lack of energy transfer in the generative region (not shown) Is 
In serious disagreement with known features of large scale turbulence 
(Lumley and Tennekes, 1970» p. 271). It remains an unresolved problem 
to devtlop a model based on a W-’H representation which will adequately 
reproduce the spectral flux convergence of the generative region. 

In summary. It Is concluded that for a realistically large Reynolds 
3 4 

number range (10 <Re<10 ) the derivative model adequately represents 
the low order moment structure and the Inertial properties of the u 
component of turbulence In the viscous subrange. For some applications 
such as those Involved with systems response at scales comparable to the 
viscous subrange, the derivative model, as described above, may be 
sufficient. But for problems Involving eddies with dimensions of the 
order of the energy containing eddies, neither the velocity nor the 
derivative models adequately represent the non<~llnear features of that 
region as manifested In a spectral transfer of energy. 

We now move to an application and evaluation of the simulation of 
realizations of turbulence by some of the linear and non-linear models 
discussed In Chapters 3.0, 4.0, and 5.0. 



110 


6.0 DEMONSTRATION AND EVALUATION OF MODELS 

The development of Chapters 3.0, 4.0, and 5.0 has concentrated on 
the calculation of the response functions from given statistical data. 
The linear filter Is derived from a predetermined spectral form whereas 
the non-llnear filters of Chapter 5.0 require additional empirical data 
In the form of third and fourth moments. It Is the purpose of this 
chapter to briefly describe the Implementation of the results derived 
In the previous chapters for the simulation of turbulence. The dis- 
cussion Is limited to linear and non-llnear models of a single velocity 
component . 

6.1 Generation of White, Gaussian Series 

Standard computer sub-routines exist for the generation of randomly 
ordered, Gaussian series. The difficulty with using these random 
number generators lies In their small but significant delvatlons from 
a Gaussian distribution for moments of order greater than 2. The 
white spectrum condition Is Improved Iteratively by several random order 
shufflings of the Input series. 

The problem of non-Gauss lanlty was overcome by generating random 
values from the cumulative probability density function of a Gaussian 
process. For a Gaussian process with mean 0 and variance 1, the 
probability, p, that a sampled value, s, will be less than x, 

p ■ Pr(s £ x) (6.1.1) 


Is given by 



p « erf(x) 


( 6 . 1 . 2 ) 


The range of p Is [0,1] for a domain of x Conversely, the 

Gaussian distributed variable x, which occurs with probability p Is 
given by 

X » erf ^(p) (6.1.3) 

The generation of N values of x from (6.1.3) was achieved for an equl- 

spaced partition of the range of probability, p. Into N Increments. The 

limit of accuracy of the routine for the Inverse error function 

(Abramowltz and Stegun, 1964) of 5 10 ^ for single precision compu- 

A 

tatlons limits the length of generated time series to about 2 10 points. 
Conversely, because the approximation to Gaussian moments of order 
greater than 2 becomes Increasingly dependent on several rare large 
deviations. It Is necessary to generate a minimum number of points to 
achieve approximate Gausslanlty In the higher moments. The minimum 
then depends on the degree of accuracy desired In the moments of the 
filtered series. The conditions on the length of the series are less 
stringent for linear simulations because of a lack of Interaction 
among moments . 

6.2 Linear Model Test 

Several methods are available for the evaltxatlon of the trans*- 


formation 
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where Is the filter and ^ and y are the white and filtered series 

respectively. The direct evaluation of (6.2.1) by either a discrete 

analog or a quadrature scheme results In an additional shaping of the 

*^1 2 

Input spectrum In addition to the filtering by |k | . Therefore, It Is 
advantageous to utilize the Fourier transform equivalent of (6.2.1), 
given by 


y(w) - urioi) ?(o)) (6.2.2) 

The evaluation of K In (6.2.2) may be made In either of several ways. 

^1 

The first method Is the direct evaluation of K from K. However, this 

method was found to be Inexact In specifying the low frequency (wave 

number) spectral content. This error In filtering the large scales 

arises from the recursion Involved In estimating the Laguerre coefficients 

(Appendix B) . A second method which Is more exact Involves by-passing 

the redunant step of computing and Its Fourier transform. The 

method of spectral factorization (Chapter 2.0 and Appendix B) results 

^1 

In an exact estimate of K directly from the given Input spectrum. 

However, In order to accommodate the fast Fourier transform techniques 

'"I 

It Is necessary to Interpolate the estimate of K , or more correctly. 

Its regularized spectral equivalent to equal Increments of o) (or f ) . 

^1 

As a result of the spectral factorization K Is tabulated at equal 
Increments of u where 

0 ) “ tan u/2 

. ' f 


(6.2.3) 
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The procedure, used to estimate K (u) was to Interpolate the regularized 
spectral factors for the u derived from the Inverse of (6.2.3) and then 
to transform the o)'-'Space using the relationship (Rlno, 1970) 

K^(co) - g(u)) (1 + expC-lu))]""^^ (6.2.4) 

where g and n are defined In Appendix B. The Interpolation scheme used 
to Interpolate g was a third-order spline function routine. The success 
of this second method Is guaranteed by the fact that the regularized 
spectral density functions are smooth, slowly varying functions of 
scale. 

✓s 

The combined error of spectral factorization to estimate g and Its 

^1 

Interpolation to estimate K was found by computations to be less than 
-3 

10 of the modulus at a glven^ scale. Therefore In estimating a 
spectrum the only discernible disparity between an Input spectrum of 
turbulence and Its simulation lies In the deviation of the white noise 
spectrum from unity or In Its statistical variability. This non- 
whiteness can be eliminated In. spectral comparisons by normalization 
of the output spectrum by the Input spectrum. Figure 6.1 presents a 
comparison of the empirical and simulated spectra of vertical velocity 

3 

for a total length of simulation equivalent to about 10 Integral 
scales. The empirical form chosen for the spectrum of w Is discussed 
In Chapter 3.0 and Appendix A. It can be seen from Figure 6.1 that the 
Kolmogorov spectral form for f » 1 has been faithfully reproduced. It 
Is concluded that within the extremes of statistical variability ex- 
pected In a single finite length of record, the spectrum of the 
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simulated turbulence normalized by the sample white spectrum Is Identical 
to the prescribed, empirical spectrum. 

This completes the demonstration of the linear filtering method, 
and we turn to an evaluation of a non-llnear model. 


6.3 Non-Linear Model Test 

The specification of a white Gaussian Input for stimulation is 
more critical for a non-llnear simulation. Deviations from a white, 
Gaussian Input process result In spurious correlations between the 
kernels of different orders. For example. In Section 5.1, If the 
input was non-white. or non-Gausslan, the moment expansions developed 
there would contain various correlation functions of the input process 
In place of the Integer coefficients. In particular, for an Input 
series, h(t), with correlation functions 


R^Ct) * h(t) h(t + t) 

R^CTi “ h(t) h(t + T^) h(t + T2) 

R^ (t^, T2» Tj) “* • • • • 


(6.3.1) 

(6.3.2) 

(6.3.3) 


the variance expansion Is given symbolically 


u^ • I • I K^K^R2 + I • I + j 

+ J • I kVr3 + J • I kVr^ + J 


I 

j K^K^Rg 


(6.3.4) 
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to third order. The lack of orthogonality (relative to a Gaussian 
p.d.f.) results In the last three terms. The non-whiteness of the Input 
process also results In multiple moment-kernel Interactions. Con- 
sequently, deviations In output statistics for a non-ideal Input are 
attributable to several sources of error simultaneously. 

Just as In the development of techniques for Implementing the 
linear model. It Is advantageous to first compute the Fourier transform 
of a simulated realization of turbulence. The expansion of the realiza- 
tion 


u(k) 


K^(k) Hj^(k) 




(p^ k-p) k-p) dp 

3 r ■ 


(6.3.5) 


+ 


# • • • 


requires only a given Input series because the are related to 
by the Gram-Schmldt orthogonallzatlon procedure. For example 


H2(k^, k^) H^(k^) H^(k2) - 6(k^ + k2> (6.3.6) 


Further, the expansion (6.3.5) reduces to 


u(k) 


y(k) + A„ 


y(p) y(k - p) dp + . . . 


+ A 2 K^(o) K^(k) + ... 


(6.3.7) 


,T. 


where 
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y(k) - ie-(k) Hj^(k) (6.3.8) 

/ 

for separable kernels. 

^e form of (6.3.7) represents a functional expansion In y with 
corrections for non-orthogonallty . The generation of y was described In 
Section 6.2. The actual simulation of u(x) Involves an Inverse Fourier 
transformation after the computation of (6.3.7). 

Some results of a sample computation for a non-Gausslan process 
with a given skewness 0.23 and a given kurtosls 3.6 are displayed In 
Figure 6.2. The method of generation of the white, Gaussian Input 
series was the same as described In Section 6.2. The procedure con- 
sisted of generating repeated, serial samples of lOil up to an arbitrary 

3 

upper limit of 10 £. The results for the skewness and kurtosls are 
not exactly those of the Input. An examination of the higher moments 
of the Input, for n ^ 5, revealed that the number of points In each 
of the Individually generated realizations of length lOi^, had a weakly 
non-Gausslan structure (C^ * 0.06, and - 14.5 rather than their 
Gaussian values of - 0 and ^ * 15). The error In the limiting 
values of skewness and kurtosls are believed to arise by spurious 
kernel-moment Interactions as discussed earlier. Two methods are 
available to eliminate such errors — either a trial and error method 
of varying the specified Input skewness and kurtosls, or a generation 

3 

of a longer, single realization, say 10 £>, so that the errors In the 
Gaussian generator are reduced. 

A second characteristic of the non-llnear test slmulatloiv Is dis- 
played In Figure 6.3. The spectrum weighted by wavenumber raised to a 




X 

Figure 6.2. Comparison of moments of simulated process for 
successive Increments to length of series with 
input values. 







1“ 



Figure 6.3. Comparison of spectra of empirically specified 

model and of simulated turbulence for a lion-linear 
representation . 
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5/3 power displays a slight systematic underestimation of the spectral 
content of the Inertial sub*-range and an overestimate at scales 
f < 1. Apparently, the spurious moment Interactions which result In 
an underestimate the total variance (Figure 6.2) vary In scale, as 
might be expected. Their accumulated effect Is to reduce the effective 
value of the response function near x ~ 0 and Increase the response at 
large scales. Further simulations utilizing single large arrays for 
Individual sample realizations can be expected to overcome this 
Implementation difficulty. In that the moment-kernel Interactions will 
be reduced. 

It Is concluded that the mechanics of simulating turbulence 
require a close scrutiny of the Gaussian and spectral properties of the 
Input process In order to achieve realistic results. Further, It Is 
recommended that computations be conducted In phase space In* order to 
utilize the accuracy of the spectral factorization procedure, with a 
transformation to real space after filtering. 
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7.0 SUMMARY AND CONCLUSIONS 

We have examined the method of simulating turbulence by filtering 
a white noise process. The development of models has been directed 
towards overcoming three areas of weakness In previous filtering 
methods. First, the calculation of filters by digital spectral factori- 
zation of enqplrlcal spectra of surface layer turbulence eliminates the 
need to use spectra with known spectral factors. Second, the develop- 
ment of linear, multi-component models follows from the methods of unl- 
component models after dlagonallzatlon of the spectral matrix. Third, 

I : : 

the non-Gausslan structure of the velocity and derivative of a 
turbulent cocqponent which Is associated with the patchy nature of 
turbulence was examined In several non-linear models. The results of 
each of these three areas of model development are reviewed and 
summarized. 

7.1 Linear, Unl-Component Models 

Because the linear filtering technique has enjoyed considerable 
success In a variety of applications, similar models were developed 
and expanded In Chapter 3.0 to Incorporate aspects of boundary layer 
similarity and general spectral forms. A comparison of several 
empirical formulations, for surface layer spectra and the kernels 
derived from them, showed that the response structure, particularly 
for large scales, was sensitive to the original empiricism. A 
definitive formulation for the energy containing sizes of the 
turbulence awaits further experimentation and theoretical consolidation. 



The Initial response and width of appreciable response were also 
examined relative to the thermal stratification. Initial response of 
velocity, normalized by surface stress, was shown to Increase In any 
stratification. Also, the range of response broadened with decreased 
stability. The simulated turbulence has a somewhat larger memory than 
the more common first order linear model — a property of Interest for 
control problems. The Initial response was shown to Increase with 
height under all stratifications, but to Increase more rapidly In un- 
stable stratifications. The response function for the vertical 
velocity component was. approximated empirically by the relationship 

K^(x;?) - exp - (x/B^)^^^ H(x) (7.1.1) 


where 

- 0.75 (1 + 0.75 |?|) 

- 0.7 (1 + 0.755 + 3.05^) 
H(x) ■ 0.5 (1 + exp - 2.5 x^^^) 


(7.1.2) 

(7.1.3) 

(7.1.4) 


over the range 
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0 < X < 10 (7.1.6) 

Further studies of the response structure as Influenced by 
viscosity Indicated a delay In the bulld-up of response to an Ivopulse. 
The delay decreased and the maximum response Increased with Increasing 
Reynolds number. The viscous model for large Reynolds number converged 
to the "Invlscld" model. The response of the turbulence In the viscous 
models resulted from the lower energy content at scales comparable to 
the viscous subrange. The memory of the simulated turbulent process 
decreased with an Increase In the Reynolds number. 

The predictive function of the turbulence Is essentially a 

single "now" value, and the mean square error of prediction grows 

rapidly with prediction distance. For example, the root-mean-square 

error of the prediction at half an Integral length Into the future Is 

estimated to be 0.8 0 , where a Is the standard deviation of the 

w w 

vertical velocity component. However, the model based on a -5/3 
spectral slope has a slightly larger memory than the usual first order 
linear model. 

A linear model of the derivative of the vertical velocity In the 
viscous subrange was examined for a later application to a non-llnear 
model applicable to the modeling of diffusion. Because the response of 
the derivative Is negative over a large range of scales, the memory 
of the derivative process Is significantly less than that of the 
corresponding velocity process. Intuitively, the reduction In memory 
Is compatible with the concept of Increased disarray In differentiated 
signals of turbulent velocity or temperature. 
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The change in derivative response with a change In the viscosity 
Indicated that the maximum response Increased and converged towards 
smaller delay with Increasing Reynolds number. The response at the 
large scales also converged to the Invlscld response, Indicating an 
Increasing Independence of the generative and viscous regions with In- 
creasing Reynolds number. Similarly, the memory of the derivative 
process decreases as the Reynolds number Is increased. 

7.2 Linear, Multi- Component Models 

The method of linear simulation was extended to the multi-component 
structure of turbulence In the atmospheric boundary layer. Specifically, 
the model was devised to simulate the lioportant cross-correlations 
between u, w, and 6. A simulation of the co-variances Is equivalent 
to reproducing the momentum and heat fluxes characteristic of the 
surface layer. 

The problem was simplified by diagonalizing the spectral matrix 
at every frequency. This was possible at scales where the matrix was 
diagonally dominant, or where the turbulence was nearly Isotropic. 

The resulting eigen-spectra were spectrally factored In the same manner 
as followed In the single component models. The response functions were 
then reoriented to the original component space by a special unitary 
transformation which preserved the minimum phase structure of the eigen- 
factors. 

The basic u-w model In neutral stratification was examined In 
terms of the corresponding single component response structure. In 
the multi-component model, the self* 7 response of u and w was found to be 
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less than the response for a single component. Although the presence 
of a vertical momentum flux decreased the self-response, the total 
response of the vertical velocity to stimuli was Increased by the cross- 
responte. At scales greater than £/10, the total response exceeded the 
sum of single conq)onent responses. Despite the Increase In larger scale 
response, there Is only a small change In the predictability of the w 
component. In addition, the prediction process remains essentially a 
single point, weighted "now" value scheme. 

The multi-component model was next extended to Include the effects 
of stratification. Consideration of the w and v components revealed 
Identical Initial responses for different variances, and different 
memories. It was concluded that the Initial response was solely a 
function of the spectral distribution of variance for scales much 
smaller than the energy-containing region. The predictability of this 
model Increases with Increasing Instability, because of the Increased 
low frequency content. The coupling between u and w also decreased 
with Increasing Instability. It was concluded that the decreased 
coupling was a manifestation of the loss of a principal axis with the 
onset of convection. The total response of w for multi- component 
stimuli decreased with thermal Instability and Increased with stable 
stratification. 

7 . 3 Non-Linear Models 

It Is well known that a linear representation of turbulence Is 
fundamentally Invalid for dynamical simulations. The difficulties In 
developing an alternative, non-llnear model are three- fold. First, the 
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degree of non-linearity varlea considerably as a function of scale, 
and second, suitable supporting experimental evidence on which to base 
even a weakly non-linear model Is scarce. Lastly, the complex Integral 
equations for the expansion of the moments must be truncated both for 
mathematical expediency and for the limited experimental Information 
about the moments. 

Computational difficulties were overcome by Invoking a simple 
closure for the relationships between kernels of different orders. 

These separability conditions are analogous to other closure schemes 
where higher order moments are expanded In terms of lower order moments. 
However, the scheme Is basically Incorrect In Its physical Interpreta- 
tion. The method of separability used In the non-linear velocity model 
assumes that the orders of response are locally Interactive whereas It is 
more realistic to expect higher order responses to be related to 
Integral properties of the lower order response. 

The advantage of the separable kernel method Is that It reduces 
the simultaneous Integral equation representation of the moments to 
simultaneous algebraic equations. These equations can then be solved 
by simple Iteration. Because the resulting expansion of the moments 
Is In terms of the linear variance, the model Is referred to as a 
quasl-llnear model. 

The results for the variance contributions of different degrees of 
non-linearity Indicated that the linear contribution was overwhelming, 
because the degree of non-Gauss Ian behavior was not large. The 
spectral content of the non-linear contributions to the variance In- 
creased with decreasing scale. The non-linear contributions erroneously 
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form an Insignificant part of the total energy even In the Inertial sub- 
range. The Initial response of the non-llnear model Is somewhat less 
than the Initial response of the corresponding linear model because of 
the adjustment of the small scale variance In terms of non-llnear 
contributions. It was concluded that a significantly larger response Is 
possible In a more realistic, non-llnear model of the Inertial subrange. 
The memory of the linear contribution to the quasl-llnear model Is 
somewhat reduced by the removal of variance at saiall scales from the 
linear part of the variance. Because the non-llnear contribution to 
the variance Increases with Instability, the memory also decreases as ^ 
Increases negatively. 

The quasl-llnear model of velocity was developed for compatibility 
with the low order moment structure of a simulated turbulent velocity 
component. However, the model was shown to be Incompatible with a 
spectral transfer of energy because of a lack of phasing between the 
kernels . 

A similar separable model was developed for the derivative process 
of the u-component of the turbulence. Estimates of the skewness and 
kurtosls as a function of the Reynolds number were based on the results 
of Wyngaard and Tennekes (1970) , The model of Pao (1965) provided an 
estimate of the spectrum near the viscous subrange. The results for the 
truncated functional expansion were found to be Inconsistent with 
measured large deviations from the Gaussian probabilistic structure. 

4 

For large Reynolds number (Re ^ 10 ), the non-llnear algebraic equations 

2 

did not converge. For Re ^ 10 , the solution of the model equations was 
physically unrealistic In that some contributions to the kurtosls were 
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3 

negative. Otherwise, for Re of the order of 10 , the separable 
derivative model Is physically realistic, and dynamically consistent 
with a non-trlvlal spectral energy transfer. The energy transfer 
function closely approximates a similar estimate based on Pao*s model 
of spectral flux convergence. 

7.4 Conclusions 

In this research, we have atten^ted to consolidate and expand the 

methods for simulating turbulence in the atmospheric boundary layer. 

The mathematical formulation of the filtering method of simulation 

has been based on empirical results of the accumulated research Into the 

spectral, spatial and probabilistic structure of turbulence. The 

simulation technique has been extended to spectra without analytical 

spectral factors, to the multi- component spectral matrix, and to the 

non-linear, non-Gausslan structure of turbulence. 

Previously, in order to produce realizations of a process similar 

to turbulence, it has been customary to truncate the functional 

expansion of turbulence at the first (linear) order, and assume a 

Gaussian probabilistic structure. In addition, the spectrum of the 

-2 

atmospheric turbulence has usually been represented as f rather than 
-5/3 

f over the range of scales of engineering interest. Although the 
resulting misrepresentation of variance may be rather Insignificant 
between these spectral models, it Is not known a priori what effect the 
differences in variance distribution will have on, say, the predictive 
structure of the turbulence for purposes of control and stability. 

These questions have been examined In terms of some aspects of large 
scale response, and of predictability. 



I 


129 

The necessary computational technique for all the model develop- 
ments was the digital spectral factorization technique. It Is 
concluded that the range of applications of simulated turbulence has 
been broadened by the numerical computation of the spectral factors of 
rather general spectral forms. 

The development of simulation models was extended to multi- component 
representations. The methodology was simplified considerably by 
dlagonallzatlon of the spectral matrix to produce multiple single 
component spectra. Each elgenspectrum Is factored as In the unl- 
component model development. The multi-component spectral factors and 
response functions are formed by rotation of the elgenstructure by the 
eigenvectors of the spectral matrix. 

It was concluded that the total response of the multi-component 
models exceeded that of corresponding unl-component models at scales 
greater than about £/2. However, the predictability of the simulated 
turbulent process was not significantly Increased In the multi- 
component models. 

The need to Incorporate non-llnear aspects In the simulation of 
turbulence, particularly In synthesizing the relatively rare, large 
scale gusts was discussed. Methods to achieve a reasonable simulation 
of the observed patchiness of atmospheric turbulence have led to the 
Introduction of the discrete gust concept. Such discrete gust models 
specify that turbulent realizations contain superimposed eddies of an 
Invariant form, such as a ramp function. The strength and frequency 
of occurrence of such eddies Is formulated for consistency with a 
Kolmogorov spectral structure and enq>lrlcal estimates of the exceedance 
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statistics of turbulence. While such methods are more realistic than 
linear, Gaussian spectral models, their dynamical basis is controversial. 
Because the invariant eddy forms, for example the ramp function, are 
expedient mathematical idealizations, there is no direct relationship 
of the parameters of the particular representation with other aspects 
of the structure of surface layer turbulence. The implications of 
such models to known dynamical and statistical structures needs to be 
studied. 

The results here are a compromise between the discrete gust methods 
and the well-known filtering methods. The minimum-phase response 
functions constitute invariant eddy forms of random strength and 
occurrence. The amplitude of each eddy is determined by the history 
of the Gaussian, white noise input realization. The kernels which are 
the discrete gust form, have been shown to evolve as a natural property 
of the spectrum of the process. The kernels also correspond to a 
deflneable physical mechanism, that is, as a response to an Impulse. 
Moreover, the response functions have been demonstrated to be funda- 
mental in describing the predictive and probabilistic structure of 
the turbulence. It is concluded that the development of discrete gust 
models is not Independent of the spectral or filtering method, and 
can be made compatible by introducing the concept of the response 
structure of the turbulence. 

The viewpoint taken here has been that the requirement for simu- 
lating large gusts can be met by systematic recourse to theoretical 
and empirical formulations of the statistical structure of turbulence 
in the atmospheric boundary layer. This approach Involves use of the 
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spectral sliollarlty theories of Kolmogorov (with appropriate extensions 
In the generative and viscous subranges), and use of published empirical 
data for the multl~component spectra against a framework of the Monln- 
Oboukhov similarity theozry. In addition the formulation of a simple 
non-llnear model of a velocity coiiq>onent has Incorporated low order 
moment data so as to be consistent with Monln-Oboukhov and Reynolds 
similarity. 

The reduction of the non-llnear simulation problem to workable 
proportions requires several mathematical simplifications. The (In- 
finite) functional expansion of the moments Is truncated to the highest 
order of available data. An arbitrary specification of the Inter- 
relationships among kernels of various orders reduces the numerical 
complexities from one of solving simultaneous Integral equations to 
simultaneous algebraic equations. Although the separable model for 
velocity In the form chosen misrepresents the dynamics of the turbulence. 
It was concluded from an examination of the separable derivative model 
that It Is possible to produce a realistic one-dlmenslonal energy 
transfer locally. It remains to develop a model based on spectral 
separability for the generative region. 

It was concluded from an examination of test realizations based 
on a linear and a non-llnear model that, within the limitations of 
generating white, Gaussian noise on a digital conq>uter and the 
implementation of numerical methods In forming convolution that the 
method is practical and forms a useful representation of the statistical 
structure of turbulence In the atmospheric boundary layer as It Is 
presently understood and described. 
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APPENDIX A EMPIRICAL RESULTS OF BOUNDARY LAYER STRUCTURE 


Some recently reported observational studies (Businger et al, 1971; 
McBean et al, 1971; Busch and Larsen, 1972) have significantly reduced 
the error variance within and between the estimation of similarity 
functions measured at different sites. For convenience, we choose the 
results of the Kansas experiment reported on by Businger et al, (1971) , 
Wyngaard and Cote, (1971a), Wyngaard et al, (1971b), Kalmal et al, (1972), 
Wyngaard and Cote, (1972) because of their extensiveness and their 
internal consistency. 

The similarity functions for shear, temperature gradient and 
dissipation of kinetic energy are given empirically 

- for unstable stratification (- 1 < C < 0) 

- (1 - 15 (A.l) 

m 

- 0.74 (1-9 (A.2) 

♦g(?) “ (1 + (A.3) 

- and for stable conditions (0 < ^ < 0.5) by 

- 1 + 4.7 C (A.4) 

♦ijCC) - 0.74 (1 + 6.4 O (A.5) 

♦g(C) - (1 + 2.5 c3/5)3/2 


(A. 6) 
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Because the 'dissipation* of tenq>erature variance was not measured 
during the Kansas experiment. It was necessary to assume 

• (A. 7) 

The second moments have been fitted as quadratic polynomials In to 
the data listed In the catalogue of the Kansas data (Izuml, 1972). 

The results are displayed In Figures A.l to A. 5. The empirical 
representations of variance and co-variance are given by the following 
empirical formulae 

■ . i 

- for unstable stratification (- 1 < C < 0) 


u^ - 6.25 (1 - 0.40 (A.8) 

v^ - 3.25 (1 - 1.0 (A.9) 

w^ - 1.35 (1 - 0.90 (A. 10) 

0^ - 3.24 (1 + 0.87 C + 0.40 C^) (A.ll) 

u6 - 3.70 (1 + 2.8 C + 2.6 zh (A.12) 


- and for stable stratification (0 < C < 0.5) 


u^ - 6.25 

(1 - 0.40 C)^ 

(A. 13) 

v^ - 3.25 

51 . ■ 

(1 - 0.14 

(A. 14) 
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- 1.35 (1 + 3.2 C - 3.0 C^) (A.15) 
0^ - 3.24 (A. 16) 
u0 - 3.70 (1 - 2.0 ? + 2.5 (A.17) 


The over~bar tilde denotes a variable scaled by or as appropriate. 

The limited range of C (~1 < C < 0.5), over which the empirical 
functions above are valid, excludes the free convection region of ex- 
treme, negative C and the extremely stable region. In the later case 
any turbulence that exists In the stable region Is considered to be 
associated with gravity waves, and its structure to be dissimilar to 
the structure In the stability range with sub*rcrltlcal Richardson num- 
ber (- 0.20) (Arya, 1972). These relationships therefore define the 
statistical structure of the turbulence In a region where thermal 
effects do not overwhelm mechanical effects. Also, estimates normalized 
by T^ for C ■ 0 necessarily represent only an average about neutral 
s tr at If lea t Ion . 

Results for the similarity structure of spectra of velocity and 
temperature within the constant flux layer have been reported by many 
Investigators (Flchtl and McVehll, 1970; McBean and Miyake, 1972; 
Panofsky and Mares, 1968; Fond et al« 1971) In addition to the Kansas 
Investigators. The normalized velocity spectra C) are usually 

represented empirically In a form 


A..(C) f 

(1 + (B i(C) f)*')^'^^ 
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Estimates of n vary between authors, the most popular being n * 1 or 
n ■ 5/3. Because there Is no clear concensus of any preferred 
eoqplrlcal form, we have chosen a slightly different form given by 


f G^^(f) 


. . 2/3 .2/3 


1 + f) 


4/3 


(A. 19) 


The implied 2/3 slope of the energy distribution for B^f « 1 disagrees 
with Kalmal's results but Is more In agreement with Busch and Larsen's 
findings. However, as It Is doubtful whether Monln-Oboukhov similarity 
holds at such large scales, the correct formulation of the spectra 
would require modifying more than the spectral slope for f « 1. The 
functions and B^^ are determined for coiq>atlblllty of (A. 19) with the 
known asyn^>totlc structure of the turbulent spectra In the Kolmogorov 
range and the variance estimates given In (A. 8) to (A. 17). The 
condition for the variance Is 


^ 2 


u 


1 


>00 

f df 

j 

0 


. 2/3 . „ -2/3 

'*’e ®li 


(3w/4) 


(A. 20) 


and, the condition for the spectra for 1 << f << *K0L’ 
the normalized Komogorov scale given by 


*K0L ■ 


Re 


3/4 . 1/4 


(A. 21) 


becomes 
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f5/3 „ . . * 2/3 . -4/3 

£ G^^<£) - 


(A, 22) 


a.. <t> 
i 


2/3 


where Is the Kolmogorov constant appropriate to the 1 variable for 
a non-dlmenslonal spectral representation. From Equation (A. 20), (A. 23) 
and (A. 24) 


A • a (“ -^)^ 4> 

11 l'‘3 ^ 

r ~2 

, . .4 "l ,3/2 . -1 

®11 % *£ 


(A. 23) 


(A. 24) 


For a value of Kolmogorov's constant, a ~ 1.5, for a three 
dimensional kinetic energy spectrum 


- (1, 4/3, 4/3) (18 a/55) for 1 - 1, 2, 3 (A. 25) 


Similarly, the empirical functions for the texoperature spectrum are 

1/3 

given by Equations (A. 25) and (A. 26), except (j) Is replaced by 

£ EX 

and the coefficient, a^, replaced by the corresponding coefficient, 
for temperature. 3^ la not as well defined empirically as a^, with 
estimates ranging from 0.4 to 0.9 (Panofsky, 1969). For consistency, 
we have selected the empirical estimate based on the Kansas date, 

- 0.8 + 0 . 1 . 

The cospectra of stress and vertical heat flux have been shown 
(Kalmal et al, 1972; McBean and Miyake, 1972) to have a ~7/3 slope In 
the Inertial subrange. These findings are In agreement with anisotropy 
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arguments (Tennekes and Lumley, 1972). However t^e horizontal heat 
flux co-spectrum has a spectral slope which varies both with height 
and stability. Kalmal, et al, report that the slope of the u6 co- 
spectrum decreases In unstable conditions from -3 near the ground to -7/3 
at the top of the surface layer, and that In stable conditions the slope 
Is -5/3. The difference between -7/3 and -3/3 Is considered Insignifi- 
cant here, and the exponent of the cospectra Is estimated everywhere to 
be -7/3. The enq>lrlcal form chosen for the cospectra Is 


5/3 


f G 


Ij 


1 + f)' 


(A. 26) 


2/3 -1/3 

The scaling factors K , which correspond to (|) and In the 

ij c X ^ 

autospectra, have been determined enqplrlcally by Kalmal to be 
- for unstable conditions 


- KjgCt) - Rig(C) - 1 

and for arable conditions 


(A. 27) 


- 1 + 7.9 c 

(A. 28) 

Kjg(C) - 1 + 6.4 C 

(A. 29) 

Kj^g(C) - 1 + 17 C 

(A. 30) 


The empirical coefficients, and , are determined from slmultane- 

jfS:' 

ous algebraic conditions Involving the variance and asymptotic structure 
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of the cospectrum. In a manner Identical to the empirical represents^ 
tion of the spectra. The asympototic constants, (corresponding 

to and 3^ in the autospectra) were estimated on the basis of the 
Kansas^ data (and the particular definition of f used in this study) to 
be 

c^^ « - 0.14 (A. 31) 


c^Q - - 0.41 (A. 32) 

Cj^Q - 0.14 (A. 33) 


Finally, the coefficients A^^ and follow from the equivalent of 
Equations (A. 25) and (A. 26) and are given by 






where the normalized variances. 



are given as 


(A. 34) 


(A. 35) 


I 


13 



(A. 36) 



(A. 37) 


i 

This completes the empirical representation of the spectra and 
cospectra — for a given scale height and stability. It is assumed 
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in what follows that the cospectra Involving v are Identically zero at 
all stabilities because of the passive role played by v In the surface 
layer (Lumley and Fanofsky, Chapter 3.0). However, In order to conq>lete 
the (complex) spectral matrix representation of the surface turbulence, 
It remains to establish the phase between the velocity and buoyant 
fluctuations. According to J. C. %ngaard (private communication) the 
phase, determined on the basis of the quadrature and cospectra of the 
Kansas data. Is Indistinguishable fropi zero. This result which differs 
from the observation of Deland and Panofsky (1957) Is adopted for 
simplicity. 

We now move to a discussion of the observations of the non-Gauss Ian 
probabilistic structure of surface layer turbulence. If we consider a 
measured turbulent velocity conceptually as the result of a summation 
of multiple random Influences exerted on the flow prior to the obser- 
vation, It Is reasonable to expect Its probability distribution to at 
least approximate a Gaussian distribution according to the Central 
Limit Theorem. Early measurements (for example, Townsend, 1947) of 
the skewness and kurtosls tended to confirm this hypothesis. Some 
estimates for the skewnesses and kurtoses of single-point measurements 
made at Kansas (Izuml, 1972) are given In Figures A. 6 to A. 11 as 
functions of C- scatter Is dlsappolntngly large for most moments. 

Apparently, the averaging time, which must Increase with the order of 
the moment In order to adequately sample the rare events (Lumley, 1970b, 
p. 73; Tennekes and ^ngaard, 1972) Is not long enough. Visually, 
there appears to be a general Increase In all the various skewnesses 
with decreasing stability. Also the kurtosls of the longitudinal 
coiq>onent, u, and the temperature are Indistinguishable from a 




Figure A. 6. Skewness of longitudinal wind fluctuations with 
stability. 
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constant » whereas the kurtosls of the vertical cbmponent decreases 
with decreasing stability. Mbnln>-Oboukhov similarity Is apparently 
only satisfied for the skewnesses of w and 6. 

A positive skewness Indicates a greater probability of negative 
fluctuations whose absolute value does not exceed a standard deviation 
and a greater probability density of positive fluctuations larger than 
a standard deviation, therefore, the w component. In near-neutral 
conditions, has more small negative and large positive excursions than 
a strictly Gaussian distribution. This arrangement Increases In un- 
stable conditions and decreases In stable conditions, possibly reversing 
Itself near C * 0.5. The temperature process Is likewise skewed. 
Presumably, In unstable condition the more probable large positive 9 
excursions are associated with the Increased probability of large 
positive w excursions. Finally, the longitudinal velocity distribution 
1b symmetric in the range, - 1 < 5 < 0.5. 

The kurtosls Is a measure of the Integrated probability density 
In the extremes of the distribution Irrespective of the sign of the 
fluctuation. For a kurtosls exceeding 3, the process would have a 
greater concentration of large, absolute occurrences of either sign, 
greater than one standard deviation, than would a Gaussian process. 

From the data displayed In Figures A. 10 and A. 11 the vertical velocity 
and buoyancy fluctuations contain more large excursions than would occur 
In a Gaussian process. Also the u component (Figure A. 9) Is associated 
with a more regular process than a Gaussian process. With Increasing 
Instability, u and 6 Indicate no pronounced change In their content of 
large fluctuations, while the w process tends to become more regular. 
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The non-Gausslan structure In terms of the skevness and kurtosls 
as a function of stability has been approximated by the following 
relationships 

- for tins table conditions (-1 < C < 0) 


M? » 0.12 a - 4.7 C - 

- 0.34 (1 - 3.2 C - 
^ and for stable conditions 

- 0.12 (1 - 2.0 C ) 

- 0.34 (1 - 2.1 C + 
while for the entire range (-1 < 



“ 3.7 (1 + 0.09 O 

> 4 - 3.4 


2.3 zh (A. 38) 

1.3 (A.39) 
(0 < C < 0.5) 

(A. 40) 

1.9 C^) (A. 41) 

C < 0.5) 

(A. 42) 
(A. 43) 
(A. 44) 
(A.45) 
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APPENDIX B DIGITAL SPECTRAL FACTORIZATION 

The spectral factorization method, discussed In Section 2.4, pro- 
vldes a technique for determining a phase characteristic to associate 
with the known modulus (or spectrum) of a complex function, so that Its 
Fourier transform will be a one-sided physically realizable filter. 

This appendix describes the Implementation of this approach computation- 
ally. Its use Is not limited to rational spectra. The use of the 
fast Fourier transform (FFT) algorithm provides a very efficient 
method of solution. Other methods such as Wiener’s predictive 
decomposition and Fejer’s solution for the roots of a polynomial 
(Robinson, 1967b) are slow and Inaccurate by comparison. 

The turbulent process represented by the empirical spectra of the 
model Is referred to as a continuous-parameter process. The necessary 
and sufficient condition for the spectrum of a continuous-parameter 
process to be factorable Is that the spectrum be absolutely continuous 
and satisfy the Paley-Wiener condition that 

df > - « (B.l) 

1 + f^ 

In other words, the spectral representation does not converge for 
£ 00 as rapidly as an exponential, exp - f^, vrhere n ^ 1, However, 

4>(f) must nevertheless approach zero In the limit of large f. 

The application of the FFT algorithm to spectral factorization 
requires that the continuous-parameter process defined over - oo < f < eo^ 
be converted to a discrete-parameter process over a finite range 
- 7T ^ u 7T, where 
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/' 

u ■ 2 tan~^ f (B.2) 

For a discrete-parameter, process with spectrum to be factorable Its 
spectrum must be regular, or alternatively, 

r 

In du > - «> (B,3) 

■'-TT 


Rlno (1970) has shown that the discrete parameter spectrum will be 

regular only If the continuous parameter spectrum approaches zero In 

2 n^l 

the limit of large f no faster than a power of (1 + f or In the 

limit 


0 < Urn 4>(f) (1 + < «> (B.4) 

f-K>C» 

where p Is some positive Integer approximating the behavior of the 
spectrum In the limit. The coordinate transformation (B.2) guarantees 
that defined by 

i|;(u) ■ (1 + cos u) (()(tan u/2) (B.5) 


will be regular on [-tt, tt] . 

We will follow the development of Section 2.4 In Implementing 

spectral factorization, but with Fourier transformations In terms of the 

1/2 

FFT algorithm. The expansion of In becomes 


, ,2TTnx ,1/2 ^ 2^rnt^ 

In W(— )] - exp(-l 


N 


(B.6) 
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Equivalently » the a coefficients for the regularized process are 


1 1 r i./2irn.il/2 ^ >2TTntv 

a - - I In [ip(-=-)] cos (-75— ) 

^ ” n-0 ” 


(B.7) 


From (2*4.20) the minimum-phase kernel transform* %, whose modulus 
Is ip» Is given by 


/ 2^^ . / 2TO\ 

g ("Y") - exp h (“Y"> 


(B.8) 


where 


h (^) - ^ -H (-1 + (-1)" ^ 


(B.9) 


The final step Is the recovery of K, where 


K K - 


(B.IO) 


from i where 


g i* " ^ (B.ii) 

The solution Is given as an approximation (Rlno* 1970) 

K(r)<r/2a~^ E L (2t) (B.12) 

n-0 ” “ 

where Is the Laguerre polynomial of order n. The coefficient Is 
derived from g by a Fourier transformation 
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K 


* ■ jo 


(B.13) 


Also, is derived from Y > ^>y iteration 


yM . 1 lt-1 

^0 /2 ^o 


(B.14) 


. 1. (^M-1 ^ M-1^ 


j > 1 


(B.15) 


^ ^ ° 


(B.16) 




j i 1 


(B.17) 


The computer program written to factor a given spectrum was 
calibrated for a spectral form 


<Kf) - (1 + 


(B.18) 


which has a known response function (for a given N) of 


^/2 ~ N -1 

K(x) ■ exp - « X X ^ 0 


(B.19) 


- 0 


< 0 


(B.20) 


The difference between the analytical and numerical solution, for various 
values of N, are displayed in Figure B.l. The error is insignificant 
up to about X * 7 and N • 4 and is insignificant to large x as N dc“ 
creases. For N • 1, the error is of the order of the single precision 



AK'W|/K'(5rt 


L63 




accuracy. Although the response function itself is relatively snall^ 

“3 ' -V ^ / 

(^10 ), for large x, say x > 5, for all N, it is of the order of the 

error in the method. As a result any simulation of turbulence at 

scales about 10 times the energy-contalning scales using the kernels 

can be expected to be in error. However » as such scales are not of 

immediate practical concern, the algorithm for spectral factorization 

is considered accurate for most purposes. 

To be factorable, a given spectrum must satisfy (B.4). The question 

arises as to the compatibility of spectra which vary as f^*^^ in the 

limit of large f, where p is not an Integer. The approach used in this 

study has been to ensure that the largest frequency, say, 

represented by the FFT is at a scale much smaller than the resolution 

desired in constructing the filter. The nearest integer approximation 

to the parameter p in (B.4) is then computed on the basis of the 

spectral energy densities at f and f /2. Succeeding conq)utatlons 

max max 

are then inaccurate only in their representation of the spectrian near 

f^^, or in the computation of K very near x - 0. Suitable convergence 
*-5/3 

for a f spectral form was found by choosing or well 

into the inertial subrange. 
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Simulation of Atmospheric Turbulence 
by Proper Orthogonal Decomposition 

John A. Dutton and Erik L. Petersen 
ABSTRACT 

A method that produces realistic simulations of atmospheric 
turbulence is developed and analyzed. The procedure makes use of a 
generalized spectral analysis, often called a proper orthogonal 
decomposition or the Karhunen-Loeve expansion. 

A set of criteria, emphasizing a realistic appearance, a cor- 
rect spectral shape and non-Gaussian statistics, is selected in 
order to evaluate the model turbulence. 

An actual turbulence record is analyzed in detail providing both 
a background for comparison and input statistics for the generalized 
spectral analysis, which in turn produces a set of orthonormal 
eigenfunctions and estimates of the distributions of the corresponding 
expansion coefficients. 

The simulation method utilizes the eigenfunction expansion pro- 
cedure to produce preliminary time histories of the three velocity 
components simultaneously, and then, as a final step, a spectral 
shaping procedure is applied. 

Two experiments are performed, providing two time histories of the 
velocity components of 30 minutes duration. This experimental tur- 
bulence is analyzed and judged to be a realistic simulation of actual 
turbulence . 

The method is unique in modeling the three velocity components 
simultaneously, and it is found that Important cross-statistical fea- 
tures are reasonably well-behaved. It is concluded that the model 
provides a practical operational atmospheric turbulence simulator. 



simulation of Atmospheric Turbulence / 

by Proper Orthogonal Decomposition 

John A. Dutton and Erik L. Petersen 
1.0 INTRODUCTION 

The demands for a realistic simulation of atmospheric turbu- 
lence have Increased over the years because of its obvious importance 
in diffusion, aeronautics, wind-loading of structures, and all 
boundary layer processes. 

Requirements to be Imposed on turbulence simulation schemes must 
produce a compromise between the accuracy with which the empirical 
statistical structure is represented and the feasibility of the 
computational scheme. A set of criteria were suggested in a previous 
NASA report by Dutton and Deaven (1971). They are, slightly modified: 

(1) The model, through variation of internal parameters, should 
be able to simulate the various intensities of turbulence in the 
atmosphere and to provide an estimate of the likelihood of occurrence 
of each time history. This flexibility makes it possible to generate 
time sequences that approach threshold (or catastrophic) intensities 
for the systems whose response is being studied and to estimate the 
probability of failure. 

(2) The model should produce time histories that exhibit the 
sequential behavior of actual turbulence. 

(3) The model should produce signals that possess the most 

notable observed statistical characteristics of actual turbulence: the 

^non-Gaussian behavior of the density function and the exceedance 
statistics and the dependence of the energy spectrum on the -5/3 
power of the wave number or frequency over a wide range. 

Standard methods that filter a white noise process so that the 
resulting spectra resemble those of turbulence fail to satisfy most of 



these criteria. Usually a linear filter is used with Gaussian white 
noise as input, and so the resulting simulated turbulence is also a 
Gaussian process, clearly contradictory to observed evidence. 

The direct use of observed turbulence, which obviously satisfies 
the last two requirements, has only limited value because the first 
requirement is not satisfied. 

The fluid motions to be modeled are described completely by the 
Navier-Stokes differential equations together with appropriate bound- 
ary and initial conditions. Despite the simplicity these equations 
possess compared to the complexity of the motions they describe, it 
is by no means feasible at the present, nor in the foreseeable future, 
to use these equations directly in operational simulation models. 

The model to be presented in this report is based on an approach 
suggested by Dutton (1968) and further elaborated by Dutton (1969) and 
Dutton and Deaven (1971). 

Information is extracted from measured turbulence by means of 
Loeve-Karhunen expansions and is carried in the model by the ortho- 
gonal functions and the statistics of the expansion coefficients. The 
method is based on the Proper Orthogonal Decomposition Theorem. 

In contrast to most models, this model generates all three 
velocity components simultaneously, and it is found that the simulated 
time histories can meet most of the requirements stated above. More- 
over, cross-statistics between components seem to be modeled satis- 
factorily to at least second order. 


2.0 REVIEW OF SOME PREVIOUS SIMULATION METHODS BASED ON THE 
PROPER ORTHOGONAL DECOMPOSITION THEOREM 

The first attempt to apply the Proper Orthogonal Decomposition 
Theorem in the study of turbulence apparently was made by Lumley 
(1965). Before that, the basic theorem, given by Loeve (1955), had 
proven useful in the study of large-scale meteorological features 
(for example, Lorenz, 1965; Kutzbach, 1967). 

The outline for the application of the method in the simulation 
of atmospheric turbulence was given by Dutton (1968) and Dutton and 
Deaven (1969). The theory was discussed in some detail and an 
attempt was made to use the method to determine if the large gusts 
as measured at Cape Kennedy had a characteristic structure. Sixty 
large gusts for each of the components u and v were extracted from a 
turbulence record. It was found that eight eigenfunctions explained 
at least 97 percent of the variance for each component, leading to the 
conclusion that the large gusts had a characteristic structure. 

This result pointed to the possibility of simulating large gusts 
by using the eigenfunctions and by sampling the expansion coefficients 
from estimates of their respective distributions. 

Before this approach can be applied in practice, it is obviously 
necessary to determine: 

(1) how to model the turbulence between the gusts, 

(2) how large a fraction of the record shall be occupied by the 
large gusts, 

(3) how much of the total variance shall be due to the large 


gusts. 
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An answer to these questions was attempted by Smith (1971) in 
the research done under NASA contract NASA-21140 for his Master of 
Science degree. Smith tested and evaluated several simulation methods 
using the three criteria described in section 1. These three criteria 
were compressed to the following tests: 

(1) Does the simulated turbulence look like turbulence records? 

(2) Does the energy spectrum fall off as the -5/3 power of 
wavenumber? 

(3) Does the probability of getting small and large values 
exceed the Gaussian probability even though the density function is 
nearly Gaussian? 

(4) Does the exceedance statistics, as represented by the 
probability of crossing a certain level per unit time, look like those 
obtained from real turbulence? 

The following description (2. 1-2. 6) of methods and their evalua- 
tion is largely extracted from Smith (1971) . Results of tests against 
the four criteria above are given in Table 1, 

2.1 Random White Noise 

Random noise as a discrete random signal that contains equal 
energy at all frequencies can be generated using the random generator 
available in almost all computer systems. It is constructed as an 
ordered set of random variables such that for each point in time, the 
random variable is selected independently from a normally distributed 




population. 



TABLE 1 - AN EVALUATION OF SOME PREVIOUS SIMULATION METHODS 


Criteria 


Generation Scheme 

Realistic Sequential 
Appearance of the 
Time-History 

A - 5/3 Slope 
in the Energy- 
Spectrum 

Non-Gauss ian 
Density Functions 

Realistic Exceedance 
Statistics 

Random white 
noise 

no 

no 

no 

no 

Shaped random 
white noise 

improved 

yes 

no 

no 

Random noise 
with gusts 

no 

Cit is too easy to 
distinguish the gust.) 

yes 

yes 

but exaggerated 

no 

Shaped random 
noise with gusts 

no 

(as above) 

yes 

yes 

yes 

but exaggerated 

Random noise with 
variable gusts 

no 

(as above) 

no 

yes 

yes 

Shaped random noise 
with variable gusts 

improved 

yes 

yes 

yes 
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2.2 Shaped Random Noise with Preserved Phase Angles 

Because random noise clearly failed criterion 2 because of the 
flat energy spectrum, an obvious step was to take random noise and 
then shape the energy spectrum to give It the correct -5/3 slope. 

This was accomplished by taking the Fourier transform of the generated 
series and changing the Fourier amplitudes according to a predescribed 
scheme. The scheme also Included an algorithm for preseirvrlng the 
Fourier phase, because these are believed to be of some Importance 
for the Intermlttency of the turbulence (see section 2,8). 

2.3 Random Noise with Gusts 

The most pronounced failure of the method above was Its Inability 
to produce the non-Gausslan nature of the probabilistic structure. 

The method proposed by Dutton using the basic structure of large gusts 
as revealed by empirical eigenfunctions was then attempted. From the 
first eight eigenfunctions, a number of gusts were constructed and 
inserted at random. into a white noise series. By trial and error it 
was found that the best result was obtained when the gusts occupied 
40 percent of the total series, and when the ratio of variance of gusts 
to variance of white noise was 19. 

2.4 Shaped Random Noise with Gusts 

The series generated In Section 2.3 was then subjected to a 
spectral shaping as described in Section 2.2. 
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2.5 Random Noise with Variable Gusts 

The same procedure as described in Section 2.3 was followed 
except that before the gusts were inserted into the series, the 
length of the gusts was made variable by randomly expanding the gusts 
by factors of one, two, or three. 

2.6 Shaped Random Noise with Variable Gusts 

The series generated above was subjected to the spectral shaping 
process. Although this model apparently turned out as acceptable, 
it does not satisfy the requirement that it should be possible to 
generate a wide range of turbulence simulations by adjusting a few 
parameters and so the tests for sensitivity to changes in the various 
parameters would be very cumbersome. Actually the parameters in the 
model are those describing the orthogonal functions representing the 
gusts, the distribution functions of expansion coefficients, the ratio 
of gusts to total record, the ratio of variance of gusts to variance 
of total record, variations in the length of gusts, the points where 
the gusts are to be Inserted, and the spectrum to be produced by the 
shaping process. 

2 . 7 Simulation Using Empirical Orthogonal Functions 

The possibility of simulating turbulence by using empirical 
orthogonal functions to represent the entire time series was in- 
vestigated by Dutton and Deaven (1971) and Smith (1971). An 
alternative approach suggested by Hirose and Kutzbach (1969) was 
applied to a sample of nine turbulence runs each of 1024 points. 
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The approach makes it possible to obtain the eigenvalues and the 
eigenfunctions by diagonalization a 9 ^ 9 matrix Instead of a 1024 ^ 
^24 matrix as required by the conventional method. Eight eigen- 
functions a.nd corresponding expansion coefficients distributions were 
extracted and then used to simulate turbulence runs. The method ap- 
parently falls because of dependence between expansion coefficients 
and dependence between and within the eigenfunctions. 

2,8 Simulation by Manipulating Fourier Phases 

An attempt was made by Spark and Dutton (1972) to assess the 
Importance of phase angle (Fourier-phases) considerations in the 
modeling of Intermittent turbulence. The conclusion of the study was 
that intermlttency appears to be dependent on some higher order 
association in Fourier space and that any mathematical model in Fourier 
space for intermlttency would extremely involved. Nevertheless, based 
on qualitative arguments it was suggested the Fourier angles might be 
used to simulate turbulence by the following procedure : 

(1) generate a random Gaussian series, 

(2) obtain its Fourier transform and the Fourier coefficients, 

(3) form the phase angles from the coefficients, 

(4) replace the original spectrum with a smoothed -5/3 spectrum 

(5) adjust a suitable number of phase angles according to a 
preassigned schedule devised by Spark and Dutton, 

( 6 ) use (5) to find the new Fourier- coefficients and back- 
transform these to obtain the simulated turbulence. 
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Because the manipulation of the angles under step 5 essentially 
creates gustlike events in the series, this model bears strong 
resemblence to the model described in section 2.6. Unfortunately, 
it also shares with it some of the disadvantages, including the 
difficulty in assessing the sensitivity to changes in the various 


parameters . 



[ 
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3.0 THE GENERALITY OF TURBULENCE PRODUCED BY MODELS 


En 5 )irical turbulence models of the type discussed in this report 
evidently will depend on information extracted from one or more actual 
turbulence records. Consider the case where a turbulence record mea- 
surement at 6 meters height at mean wind speed of 6m/sec has been used in 
the creation of the model, and the model turns out to work satisfactorily 
in the sense of simulating turbulence resembling actual turbulence 
occurring under the above conditions. Will it then also be possible to 
use the model to simulate turbulence as found in clear air turbulence 
with, say, a length scale of 600 meters and a mean wlndspeed u = 30m/ sec? 
If the answer is no, the usefulness of the model will certainly be 
severely restricted. 

Fortunately, measurements in the troposphere and the stratosphere 
show that turbulence possesses in the inertial subrange a high degree 
of self-similarity^. Because amplitudes are related to the 1/3 power of 
the wavelength, if the wavelength in a turbulence record obtained in 
the inertial subrange is expanded by a factor 3» then an amplification 

1/3 

of the velocity amplitudes by 3 would yield a record similar to 
turbulence (Dutton and Deaven, 1969). 


^In the sense of B. Mandelbrodt ’ s (1965) concept of self-similarity : 

A process with random variable X is self-similar if the variable 
obtained by the magnification of the wavelength by h can be represented 
as a suitable magnification of the amplitude of X so that both X and 
Xj^ have the same probabilistic structure. 
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From above, for the first example, we have: a length scale of 

L = 6m and a velocity scale of u = 6m/ sec, yielding a time scale 
T = L/u = 1 sec. 

For the second example: u = 30m/sec and L = 600m, T = 20 sec. 

Thus we have to expand the temporal scale by a factor of 20, hence 

1/3 

the velocity amplitudes have to be magnified by a factor 20 

The "1/ 3 power law" is only valid in the inertial subrange where 
the spectra fall off as -5/3 in a log-log plot. In general, the 
concept of self-similarity is only useful for processes in which the 
spectra exhibit power law behavior. Unfortunately, most turbulence 
spectra in the form kS(k) have a rather flat maximum from which the 
spectra fall off towards both higher and lower numbers. 

However, for many engineering applications, it is usually the 
energy content in the inertial subrange that is of most Importance, 
and so if turbulence can be simulated in this range by an appropriate 
model, then time and velocity scales can be changed as required by 
particular applications. 

Another strategy would be to assume the spectra exhibit different 

power laws, on opposite sides of a particular wave number, K„.„, an 

MAa 

assumption which finds some support from observed turbulence. Then 
we would have 

kS(k) = Ai^ S(k) = Ak^ ^ K < 

for 

^ S(k) = k 


kS(k) = Bk 
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Denoting the temporal multiplication factor by 3 and the velocity 
factor by h, we have from Dutton and Deaven 



1/3 

where it follows that Y + 1 = 5/3 gives h = 3 ^ • 

In this approach, we would obtain the Fourier transform of the 
simulated sequence and then all amplitudes on the side K < 
would be adjusted by the factor appropriate for those wavenumbers 
while all amplitudes on the other side would be adjusted by the other 
factor. These adjusted amplitudes would be used in a back- transform 
to produce a series appropriate for the expanded time scale. 
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4.0 THE THEORY AND THE MODEL 


The experience gained through the various studies mentioned in 
the previous sections emphasizes the difficulties in generating the 
correct time sequence of the empirical turbulence, even if reasonably 
chosen statistics seem to be modeled rather well. But the experience 
also showed that one way to success could be to emphasize the creation 
of gust-like events in the generated series, thus simulating the so- 
called "surprise” of real turbulence. The current model is based on 
such aft approach, and the philosphy behind it is given a fuller treat- 
ment in the last section. For the present, it is sufficient to note 
that in order to find the structure of the gusts it is necessary to 
select a certain type of analysis and perform it on records of actual 
turbulence , 

Imagine the turbulence records to be composed of intervals of 
"passive turbulence" and of intervals of "active turbulence", where by 
active turbulence intervals we understand sequences of the records in 
which it is observed that a gust prevails for some specific length of 
time, T. Hypothesize further that the active turbulence is composed 
of a quasi- deterministic gust structure to which is added passive 
turbulence. We can then confine our attention to the active intervals 
because, if we can find the gust structure, and if we can find a way 
to represent the passive turbulence, we know from an analysis of the 
turbulence records how the passive and the active Intervals are 
distributed in the records , and a model proposes itself : generate 

series of passive turbulence and add to it the gust structure in such 
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a way that the passive and the active Intervals become distributed as 
In the actual record. 

It Is shown In Section 8 that imposing certain prin- 
ciples in order to find the gust structure and an economical 
representation of the passive turbulence lead to the Fredholm integral 
equation 


I R(s,t) dt = W(s) 

0 


where 


R(s,t) = E{f(s)f(t)} 

and f (t) Is the nth interval (t < t < t + T) of active turbulence 
n n— ~n 

taken out of the records and then redefined over the interval (0,T). 
The expectation operation E{ } is performed on the ensemble of 
f^(t) , n = 1,2, . . . 

The (|)'s are orthonormal eigenfunctions of the correlation matrix 
R(s,t) and under the assun 5 )tions stated above, reveals the gust 
structure and 4 > 2 » ••• provide us with an optimal expansion of the 

passive turbulence over intervals of length T. 

We have for an interval of active turbulence 

OO 

(4.1) g^(t) = Z a^\(t) 

k=l 


and for an interval of passive turbulence 
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o9 

(4-2) ^ ^ 


The probability density functions of the expansion coefficients 


“k 


are estimated from the f *s 

n 


n 

“k “ 


fn^t) 4>k(t) dt 


hence the in (4.1) and (4.2) are sampled from their respective 

distributions . 

Representing the gust structure by a function that is orthogonal 
to all the functions used in the expansion of the passive turbulence 
obviously requires that the two processes, the gust and the passive 
turbulence, be orthogonal, a requirement we only can expect to be met 
approximately. Then, because of a possible non-orthogonality between 
gust and passive turbulence, we would expect (j)^ to give most of the 
gust structure plus a little of the passive turbulence, ^2 some of the 
gust and more of the passive turbulence, and so on. To account for 
this in a generating scheme, one approach would be to use all the (|)*s 
to construct all sequences of the turbulence, but to diminish the 
amplitudes of the first eigenfunctions approximately in the intervals 
with passive turbulence, for example by transforming the probability 
density functions for the expansion coefficients. 


4.1 The Model 

The discussion above proposes that the model be established in 
two parts: an analysis scheme and a generating scheme. Thus, the first 

scheme describes how to obtain the eigenfunctions and the second scheme 


how to use them. 
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The analysis scheme 

The steps In carrying out the conqjutation are: 

A1 Select one or more observed turbulence records. 

A2 Select a characteristic feature in the records believed to be 
of importance for the Intermittency of the turbulence and select the 
corresponding time interval, T. (Such a fea-iure may be a large gust.) 

A3 Extract as many as possible time intervals containing the 
feature in order to construct a representative ensemble f^(t). Estimate 
the probability density function for the time interval between the 
occurence of the feature. 

A4 Subject the ensemble to a Proper Orthogonal Decomposition to 
obtain eigenfunctions and expansion coefficients and estimate the 
appropriate probability density functions of the coefficients. 

The generating scheme 

Because the eigenfunctions are defined on a fixed Interval, it is 
necessary to generate the turbulence in multiples of this interval. 

From A3 we know the distribution of the time intervals between the 
events, and we can then pick values from such a distribution to determine 
where in the generated series the special events shall occur. From the 
appearance of the eigenfunctions, one can estimate how much active 
turbulence and how much passive turbulence the first few eigenfunctions 
explain. According to this estimate the weight of these eigenfunctions 
is diminished in the passive intervals. 

The generating scheme then becomes : 

G1 Generate a set of random numbers, using the probability 
density functions from A3, to establish the sequence of intervals of 
active and passive turbulence. 
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G2 Generate the active Intervals by sanqjling the expansion 
coefficients from their probability density functions (from A4) and 
multiply the respective expansion coefficients and eigenfunctions and 
add the functions together. 

G3 Generate the passive intervals as above, but make appropriate 
transformations of- the first few probability density functions, 
and if necessary, 

G4 Obtain the Fourier-Transform, let the phases be unchanged, but 

1. Change the amplitudes such that the energy spectrum is 
proportional to the -5/3 power of the frequency over a certain range 
(see Section 7). 

2. Obtain the appropriate time and length scales by multiplying 
the amplitudes by a constant which has been determined in accordance 
with Section 3. 


Then back transform to the time domain. 
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5.0 CONSTRUCTION OF THE MODEL 


The ensemble f (t) was selected from a turbulence record that 
will be described In the next section. It is composed of 0.10 second 
block-averaged values of the u,v, and w components; the total length 
of the record is 50 minutes. 

The selection of the ensemble fjj(t) was then done heuristically 
in the following way: 

1. T was chosen to be 5 sec (=50 data points) 

2. The record was divided into 100 segments, each of the length 
30 sec (=300 data points) 

3. Inside each segment, the maximum value of w was found and a 

5 second interval of data centered around this value was picked for all 
the components 

4. The mth ensenble function f , was then obtained from the mth 

segment for 1 ^ 100, and the function was constructed by patching 

together the u,v, and w components (150 data points) sequentially. 

The occurence statistics were then very simple, and a gust interval 
of 5 seconds duration was placed at random Inside each sequence of 30 
seconds length. 

The numerical approximation of (4.3) becomes 
(5.1) ZR(s,t)(f)^(t) = 

which is the usual eigenvector equation used in principal component 
analysis. 


We have 
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R(s,t) = - Z F^(s) F (t) 
.m mm in 


where T denotes transpose and F^ is a 1 x M and 

The matrix F has the following structure 
m 


T 

F a M X 1 matrix, 
m 


F 

m 


(t) 




U (50), V„(l), V_,(2),...,V^(50), 
m m m m 

W^(50) ] 
m 


which gives R(s,t) the structure 


UU UV UW 
VU W VW 

wu wv ww 


where, for example, UV is the 50 x 50 uv correlation matrix. 

A more straightforward way to perform the analysis would be to 
create three ensembles, one for each of the components, u, v, and w, 
and then get three sets of eigenfunctions by solving (5.1) for 
R(s,t) = UW, W, and WW respectively. But then we would have disregarded 
all cross statistical information, and it would be difficult to intro- 
duce it into the analysis later. 

The way the eigenfunctions are constructed in this analysis by 
patching together u, w, and w enables us to use all the second order 
cross statistical information available to construct the eigenfunctions. 
The function 4>^(t) would then give the most likely simultaneous occurence 
of u, V, and w during a gust in w. 

The patching of ensemble functions has been used by Jaspersen 
(1971) to analysis vertical profiles of meteorological variables. 
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5.1 'Hie Eigen valtte Spectrum and the Eigenfunctions 

A diagonallzation was performed on the 150 x 150 matrix, 

R(s,t), giving the eigenvalue spectrum shown in Table 2 and the 150 
eigenfunctions of which the first 14 are shown in Figure (5.1) and the 
first 20 are listed in Table 3. 

The first two eigenfunctions seem to explain the average values 
of u, w, and w during the "gust", and it is not surprising that this 
is the most highly cross-correlated feature. We could have prepared 
the ensemble from which the eigenfunctions were calculated in such a 
way that the mean of u, v, and w was zero in each ensemble function, 
or we could have removed the ensemble mean from each ensemble function. 

In both cases we would then have to carry some additional statistical 
information in the model which can be carried by the eigenfunctions 
and the expansion coefficient distributions themselves. 

The third eigenfunction is mostly devoted to the peak in w, and 
to a lesser degree so are the fourth and fifth eigenfunctions. 

Fig. (5.1) reveals how the eigenfunctions, as the number increases, 
tend to explain features on a smaller scale. 

Table 2 shows that 72% of the total variance in the ensemble is 
explained by the first 5 eigenfunctions, 82% by the first 10, and 91% 
by the first 20. In the final generation scheme, the first 20 eigen- 
functions were applied. 

Because only 100 ensemble functions were used to create the 150 x 150 
matrix R(s,t), the actual order is 99, and only 99 non-zero eigenvalues 
can be calculated. This creates some Interdependency between the 
eigenfunctions. However, the appearance of some degree of dependence 
does not seriously compromise the method. 




186 



187 


TABLE 2 

EXPANSION STATISTICS FOR THE FIRST 20 EIGENFUNCTIONS 
Elgenfet. Expansion coe'ff. Eigenvalues 

Nr. Mean Stand, dev. % Acciim. % 


1 

-1.72 

6.63 

47.0 

27.0 

27.0 

2 

1.39 

6.28 

41.4 

23.8 

50.7 

3 

3.39 

2.92 

20.0 

11.5 

62.2 

A 

0.96 

3.15 

10.8 

6.2 

68.4 

5 

0.83 

2.29 

5.95 

3.4 

71.8 

6 

-0.22 

2.21 

4.92 

2.8 

74.7 

7 

0.17 

2.00 

3.98 

2.3 

76.9 

8 

-0.06 

1.91 

3.67 

2.1 

79.0 

9 

-0.09 

1.70 

2.89 

1.7 

80.7 

10 

0.21 

1.63 

2.69 

1.5 

82.2 

11 

0 


2.26 

1.3 

83.5 

12 



2.20 

1.3 

84.8 

13 



2.03 

1.2 

86.0 

14 



1.68 

1.0 

86.9 

15 



1.56 

0.9 

87.8 

16 



1.46 

0.8 

88.7 

17 



1.26 

0.7 

89.4 

18 



1.19 

0.7 

90.1 

19 



1.08 

0.6 

90.7 

20 



1.04 

0.6 

91.3 
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TABLE 3 - LISTING OF THE FIRST 20 EIGENFUNCTIONS 
ElSKNFUNCTlfJNS 1 " 10, POINTS lOJ -> 150 I TH£ N COHPONENT 

ElGENFiJNCTION 


1 

2 

3 4 

5 

0,04206 

Q.043B9 

-0,02797 -0,08581 

-0,00463 

0,04081 

0,04l6u 

-0,02414 -0.12440 

-0.04763 

0,04763 

0.03728 

-0,02649 -0,12406 

-0,08169 

0,03543 

0,03181 

-0,01783 -0,14562 

-0,13505 

0,03569 

0*02959 

-0,02074 -0,14946 

-0,16944 

0,03065 

0,02801 

0,00547 -0,13325 

-0.14963 

0,03531 

0,02330 

0,00501 -0,14107 

-0,14873 

0,03729 

0,01678 

-0.00723 -0,15541 

-0,11165 

0,03315 

0.01598 

0,00937 -0,20086 

-0,10016 

O,o3b90 

0,01299 

0.01624 -0,18966 

-0,08036 

0«o3226 

0,02259 

-0,00139 -0.14884 

-0,09328 

0,04470 

d, 02147 

0.02805 -0,11219 

-0,06660 

0,04119 

Q,d2d3b 

0,01994 -0,13702 

-0,06330 

0,05724 

0,02577 

0,01675 -0,13657 

-0,06202 

0,04020 

0,02001 

0,00250 -0,14097 

-0.03094 

0,04767 

0,01477 

-0.00S40 -0,13932 

o.oiloo 

0,04190 

0,02738 

0.01973 -0, 1802b 

0,02416 

0,03003 

0,d3dSl 

0.04513 -0,15510 

0.04619 

0,03914 

0,0397b 

0,05669 -0,12474 

O.OlBBi 

0,02324 

0,03649 

0,06973 -0,11510 

0,04623 

0,02599 

0,04075 

0.11781 -0,12141 

0,01166 

0,03628 

0,04732 

0,16444 -0,12159 

-0.00252 

0,04036 

0,05686 

0,18183 -0,06258 

0,01485 

0,04408 

0,05902 

0,24642 -0.03696 

0.0460/ 

0,06744 

0,07061 

0,33675 0,06246 

0,13729 

0,09402 

0,09657 

0,50081 0,17045 

0,30385 

0,07177 

0, 0898 u 

0,31521 0,12002 

0,12825 

0,05552 

0,07756 

0,16737 0,03343 

-0.02116 

0,06373 

0,06l69 

0,14509 -0,02709 

-0,03690 

0,04609 

0,0339a 

0,14624 -0,02670 

-0,07162 

0,03866 

0.04018 

0,14697 -0,02068 

-0.01957 

0,04100 

0,04274 

0,11854 -0,06016 

-0.01204 

0,04654 

0,03722 

0,13155 -0,05263 

-0.04680 

0,05133 

0,05453 

0,09696 -0,08474 

-0.12999 

0,0S076 

0,04241 

0,08998 -0,00854 

-0. 1 1681 

0,06192 

0.03346 

0,09272 -0,09174 

-0,15449 

0,06152 

0,02018 

0.06683 -0,06410 

-0.13872 

0,05577 

0.03073 

0,08516 -0,03320 

-0.12680 

0,06626 

0,06401 

0,11209 -0,07Sll 

-0,10986 

0,06197 

0,04246 

0,10604 -0,04339 

-0,13220 

0,04793 

0,04/78 

0.0911S -0,02115 

-0,09404 

0,03735 

0,03351 

0,06541 0,03253 

-0.11768 

0.03416 

0.04054 

0,06947 0,02845 

-0,13926 

0,03275 

0.03H50 

0,04315 0,00512 

-0.13950 

0.04333 

0.02138 

0,06380 0,04069 

-0,16254 

0,05680 

0,02753 

0,07046 0,05918 

-0,17244 

0,08358 

0.02323 

0.05452 0,05223 

-0,21742 

0,04834 

Q,0262o 

0,05092 0,02523 

-0,18563 

0,03762 

0,01868 

0,04071 0.0338b 

-0.17822 

0,02li4 

0,02995 

0,02220 0,01798 

-0.14597 

0,07741 

0,07562 

0.08443 -0,05637 

-0,09790 


b 7 8 9 10 

wO, 18030 0,19475 sO.lSOOb 0,04910 0,00224 

■^0,U37b 0,20585 17187 0,11482 **0,0B762 

-?0, 11044 0,22616 -^0. 20451 0,07148 -0,07436 

-?0, 12440 0,24876 ^^0, 16648 0,04803 «0, 04762 

-0,09793 0.20090 ^0,18632 -0,01538 -0,04336 

-0,03332 0.19717 ^0,15749 0,00247 -pO, 01444 

-0, 03103 0.15961 ^0, 08936 -0,01661 0,00850 

p»0. 06508 0,10041 -0.00897 -0,09484 0,13964 

^0,12310 0,05649 f?0, 03765 •*0,06899 0,15014 

-0,16348 0,05509 0,06043 -pO, 10717 0.13249 

P‘0, 12046 -0,00760 0,04969 -0,l8441 0,09698 

-0,11473 -0,05921 0,03403 -0,l9l78 0.12276 

-0,09598 -0,14827 0,00654 -0,14193 ojl22o7 

^0,07455 -0,20965 0,02609 -0,06035 0,09553 

-0,12676 -0,20969 -0,01728 -0, 02754 0,07973 

-0,10488 -0,17647 -0.02024 0.03502 0,07443 

-0,11262 -0,14657 -0,03500 0,10122 0,06765 

-0,10466 -0,13660 -0,06051 0,08827 0,10674 

-0,08477 -0,22298 -0,04703 0,13343 0,049o3 

-0,02557 -0,20068 -0,06677 0,14174 0,00397 

-0,00022 -0,i75S<l -0,12610 0,20154 0,00783 

0,00066 -0,16866 -0,17291 0,20463 0,00760 

0,04940 -0,16356 -0,14064 0,15665 0,05365 

0.05312 6,00568 -0,11573 0,16844 0,07625 

-0,01027 0,10859 -0,08119 0,03728 0,09279 

-0, 10615 0,17679 -0.03966 -0.05826 0,2lS27 

-0,06639 0, 11462 0,08750 -0,06038 0,09310 

-0,08492 0,00764 0,18096 -0, 12366 -0,03007 

-0,06446 0,00823 0,1/638 -0,16599 -0,06999 

-0,03872 0.0G516 0,14867 -0,13519 -0,00975 

-0,02856 0,04137 0,13508 -0,073l9 -0,U692 

0,03643 -0,03757 0,05377 -0,07016 -0,12204 

0,01168 -0,07781 0,03032 -0,09439 -0,ll369 

0,10268 -0,04540 0,03762 -0,09ai6 -0,09309 

0,14014 -0,07145 0,03059 0,00367 -0,04726 

0,19044 -0,02181 -0,01298 0,04520 -0, 00838 

0,2573? -0,04249 0,00599 0,06692 0,00732 

0,26244 0,02164 0,07U6 0,08718 -0,il679 

0,18101 -0,00315 0,06551 0,12372 -0,09897 

0,13255 0,00482 0,10630 0,ll302 -0,09997 

0,14719 0,02228 0,07108 0,03773 -0,l06l4 

Q,12BoB 0,05409 0,02566 0,09l92 -0,05138 

0,034/9 0,05272 0,0l0l0 0,07256 -0, 04703 

0,07750 0,02081 0,02411 0, 10115 -0,06141 

0,03805 0,06536 0,00303 0,02l26 -0,03373 

-0,03749 0.08191 0,01692 -0,03258 -0,06743 

-0,02139 -0,00500 0,01750 -0, 09941 -0,05372 

0,03201 -0,04456 0.04827 -0,12165 -0,03964 

0,01389 -0,09210 0.05976 -0,08513 -0,00726 

0,05093 -0.06935 0,03040 -0,q9305 -0,01299 

-0,05196 -O,i002l -0,08565 -0,08500 0,03746 
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TABLE 3 - LISTING OF THE FIRST 20 EIGENFUNCTIONS (continued) 
EI&ENFUNCTIONS n ^ so , points i » so i the u component 


EIBENFUNCTION 


11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

Q«03746 

0. 01031 

-0,05298 

-0,06590 

-0.04472 

0,05953 

•0.10501 

•0,03323 

•0,14053 

0,02210 

**0,0043S 

-0,00463 

-0,06376 

-0,10204 

-0,04699 

0.02974 

•0,081i2 

•6,05525 

•0,11389 

0,01297 

n,oUoa 

0,01883 

•0,07415 

-0.08231 

-0,08090 

0,03044 

-0,08794 

•0,06041 

•0,06453 

0,02259 

-0,06439 

0, 02611 

• 0,07 064 

-0,06663 

•0,08710 

0,04397 

-0,09231 

•0,06306 

•0,04795 

0,050<»2 

•0,02430 

-0,01051 

-0,07468 

-0.04378 

-0,04279 

0,04248 

•0,01053 

•0,08716 

•0,04302 

0,08855 

-0, 023S2 

-0,04473 

-0,O5fi«»2 

-0,02928 

0,00446 

0,02B15‘ 

0#00858 

-0,05748 

•0,04433 

0,65305 

6,01000 

-0,03701 

-0,06614 

•0,01594 

0.01993 

o,02is:i 

0,00102 

•0,05226 

•0, 02237 

6,65226 

0,00589 

-0 , 06322 

-0,06208 

-0,00767 

0,02940 

0,00537 

0.02417 

•0,04158 

0,01296 

0,05817 

0,00273 

•0,05462 

-0,06944 

0,02714 

0.09904 

0,02258 

0.02434 

0,05206 

•0,01333 

6,02258 

-0,01036 

-0,00059 

-0.07365 

0,01670 

0,09098 

•0^01417 

•0,02239 

0,05112 

0,04225 

-0,00176 

-0,01453 

-0,02384 

-0,06264 

0.02775 

0,05915 

•0,02352 

0,00657 

•0,00718 

0,04510 

6,00344 

-0, 01516 

-0,00408 

-0,05142 

0,02605 

0.02845 

•0,02854 

0.03126 

-0,01733 

0,02830 

0,00775 

•0,02275 

0.01852 

-0,01772 

0,05782 

0.04657 

-0,04481 

0,05566 

0,00653 

0,03840 

•0,64218 

-0,02563 

0.03387 

-0.00428 

0,06575' 

0.05445 

•0,05625 

0,04018 

0,03877 

0,03764 

•0,00143 

-0,03542 

0,02865 

0.03225 

0,03546 

0.05278 

-0,04428 

0,01563 

6,04369 

0,07664 

0,01103 

-0,03143 

0,01493 

0.05736 

0,04-013 

0.01975 

•0,08385 

•0.03397 

0,06571 

0,07343 

0,02603 

-0 , 05006 

0,0213b 

0.04647 

0,01025 

-0,0i525 

•0,11852 

-0,02082 

0,05979 

0,01572 

•0,00714 

-0,02438 

0,05134 

0.03299 

0.017B8 

-0.02197 

-0,06826 

0.00668 

0,04149 

0,03096 

•0,02038 

-0,03267 

0,06409 

0,07072 

0,03422 

-0.04964 

-0,06621 

-0,02444 

0,06094 

0,03775 

•0.0510$ 

-0,07707 

D ■ O34&4 

0.07472 

0,05084 

-0,08811 

-0,07000 

-0,02089 

0,08591 

0,07516 

-0,07695 

-0,07113 

0,01988 

0,07540 

0,03792 

-0.11230 

•0,05088 

-0.03620 

0,13771 

0.04036 

•0.09611 

-0,08689 

0,01724 

0,07527 

-0,01427 

-0.12246 

0,00213 

•0,04078 

0,18394 

0,02035 

-0,69256 

-0,11777 

-0,03238 

0,08256 

-0,06252 

-0,09488 

0,03204 

0.00772 

0.14560 

•0,00612 

-0,67492 

-0, J6274 

-0,0647b 

0. 11526 

-0,05782 

-0,03959 

0.06154 

0,02929 

0,05599 

•0,00223 

-0,06560 

-0,10425 

-0,10085 

0.11144 

•0,08713 

-0,03549 

0,06415 

0,01589 

0,04667 

•0,04937 

•6,03960 

-0,05350 

-0,09233 

0,06832 

-0,11699 

0,00683 

0,07781 

0,00619 

0,00985 

-0,00540 

-0.68689 

-0,04512 

-0,10315 

0,03959 

-0,14175 

0,06164 

0,06680 

0.00786 

-6.02983 

•0,07226 

•0.02046 

-0,03411 

-0,10866 

0.03713 

-0, iS79B 

0,11695 

0,10293 

0,02720 

-0.02593 

•0,04906 

0,01795 

-0,00402 

-0,14975 

-0. 00228 

-0.15052 

0.12611 

0,0984b 

0,00939 

-0^03671 

-0,03620 

0,08902 

0.06266 

-0,14876 

0,00442 

-0,10668 

0.U33 3 

0,03958 

0,00067 

•0,00179 

-0,06041 

0,08952 

0,07288 

-0.12864 

0,01179 

-0,05907 

0.08901 

-0,05997 

0,03436 

-0,11034 

•0,05436 

0.14692 

0,09763 

— 0 « 05765 

-0.OJ844 

-0.04035 

0,103)8 

•0, 10504 

0.069.55 

-0,13006 

•0,03774 

0,15546 

0, 13q60 

-0,01504 

0,01720 

•0.00795 

0.07171 

-0,11520 

0,06601 

•0.13076 

-0, 05661 

0,08077 

0.16960 

0 , 0 U 6 1 

0.02676 

0.03894 

0,02729 

-0,10533 

0,00688 

• 0.U636 

-0,07000 

0,06944 

0,14875 

0 , 0 5 6 4 1> 

0.01377 

0.07752 

-0,01376 

-0,08651 

-0.01630 

•0,02632 

•0,06169 

0,09637 

0.14220 

0,05965 

0,03454 

0,04767 

-0,03867 

• 0,095.38 

-0,04348 

•0.02924 

0,01656 

0,06503 

0,09532 

0,06332 

0.02222 

0,03873 

•0,04165 

-0,07756 

• 0 , Cl 0 2 / 4 

0,01660 

0,02212 

0,02997 

0.o7034 

0,04892 

0,01301 

0,03648 

-0,00951 

-0,02080 

-0,02424 

•0^00526 

0,05537 

0,00202 

0,07551 

0,05506 

0.03979 

0,04316 

-0,04084 

0,02151 

•0,02555 

0,01999 

0,07665 

0,01552 

0,09736 

0,06976 

0.03935 

0.03462 

-0,06004 

0,06292 

•0,06500 

-0,03600 

0,00675 

0,00722 

0,10863 

0,08189 

0,06507 

0,07469 

-0,06656 

0,06292 

•0,07352 

-0, 01037 

0,07746 

•0,01015 

0,10106 

0,08123 

0,05449 

0,07736 

-0,04892 

0,10642 

-0,00567 

0,04660 

0,0?916 

0,00073 

0,05957 

0,05774 

0,01795 

0,0289b 

-0,03011 

0,06025 

0,02194 

0,06576 

0,05341 

0,00546 

0,02741 

0 • 0 1 25q 

0.00500 

0,07176 

-0.02286 

0,02603 

0»030i6 

0,03165 

-0,02326 

0,01971 

-6,00341 

-0,03623 

-0.05288 

0,05211 

-0,00265 

0,02052 

0,00711 

0,06050 

•0,16405 

-0,00926 

-0,t)3SB7 

* 0 , 6 V 7 5 7 

-0,07695 

0.05815 

0,01550 

0,02936 

0,03385 

0,05177 

•0,07633 

-0,01550 

-0,0^495 

•0.06603 

-0.04863 

0,04420 

0,03069 

0,06145 

0,06429 

-0.00750 

•0,04485 

-0,03215 

-0,09049 

-0.07368 

-0,06369 

0,06174 

0,03204 

0,02548 

0,036b] 

-0,00175 

0,00047 

•0,15006 

-0,06301 

-0.04907 

-0,05933 

0,07893 

0,03448 

•0,01214 

0,05033 

•0,01416 

-0,00070 

•0,21790 

-0,o763B 

•0,05314 

-0,05424 

0,06402 

0,00924 

-0,01971 

0.05578 

0,02764 

•0,02089 

•0,23445 
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TABLE 3 - LISTING OF THE FIRST 20 EIGENFUNCTIONS (continued) 
EIUfcNFUNCriOMS 1 - 10, POINTS Hi - l«0 t THE V COMPONENT 

EIGENFUNCTION 


1 2 3 0 y 

-0,07000 0,0U-«3 -0,04009 0.07408 

0,07004 -0,10-iUu O.OISQS -0,04789 0,0054^ 

0.0t»8i5 -U,i(ib!>? 0.02303 -0,02450 0,08692 

0,06600 -Q.lliSl 0.01634 -0,0ly^‘0 0.08428 

0,07641 -0.11062 0,01656 -0,01633 0.07520 

0.07682 -0,10009 O.OPObS -0,02025 0.06798 

0,07795 -0,10954 0.01203 -0,01492 0,0/480 

0,07716 -0,l069j. 0,01070 -0,01 12() 0,05666 

0,07765 -0,10757 0.00913 0,00023 0,04493 

0,06298 -0,10564 0,00325 -0,01953 0,05065 

0,06425 -0.10554 0,00306 -0,02990 0.04939 

0,08345 -0,i079fe -0,0C2l6 -0.02764 0,03505 

0,1)8550 -0,10946 0,00666 -0,01874 0.02710 

0,06579 -0,11007 0.00405 -0,01544 0,02909 

0,00673 -0,11034 -0,00239 -0,02428 0,03966 

0,08571 -O.UOOB 0,00013 -0,03564 0.04326 

0.08485 -0,11837 -0,00086 -0,04288 0,0320/ 

0,08940 -0.114/7 -O,00303 -0,052l6 0,03030 

0,09235 -0,U&52 -0,00256 -0.05067 0,007«4 

0,08804 -0,11904 0,03367 -0.04745 0,00266 

0,08938 -0.12079 0,01648 -Q,04}6? -0.0024ft 

0,09111 -0,12858 0,01513 -0,01969 -0.02460 

0,08351 -0,12629 0,00720 -0.01335 -0,01254 

0.06137 -0,12155 0.00648 -0,01850 -0.G260H 

0.08795 -0,12949 0,01765 -0,01540 -0,05473 

0,06657 -0.12619 -0.01019 -0,01407 -0.0B232 

0,08534 -0.U727 -0,01260 -0.015/7 -0,04513 

0,08384 f0.it47o 0,00359 -0,01735 -0,05715 

0,07834 -0, 1 1 919 0,00064 -0,00070 -0,06458 

0,07788 -0,12765 -0,00005 0,01170 -0.04244 

0,0763ft -0.1227b 0,00139 0,02177 -0,03590 

0,0609ft -0,11837 -0,00139 0.03336 -0,03173 

0,06001 -fl.il720 -0,00167 0,0236ft -0,00606 

0,0/265 -(-,ilB29 0,0(JB97 0,01 ?0l -0,00410 

0,07269 -O.iliib -0,00115 0,01335 0, 0124ft 

0,0/069 -0. 11340 -0,00292 0,01564 0,010// 

0,07207 -0,1 1604 -0.00656 0.01969 0,«U5'3 

0,07166 -0,12392 -0,00635 0,0159ft 0,01367 

0,07301 -0,12412 -0,00246 0,00857 0.0243/ 

0,07354 -O.U47ii 0,00667 0, 00806 0,02366 

0,06314 -0.40596 0,02623 0. 0002ft 0.OU13 

0,07210 -0,10553 0,02955 -0,00337 -0,01222 

0,07076 -0,10069 0,02562 -0,00320 -0,02487 

0.07930 -0,10506 0,02620 0, 00620 -0,02501 

0.07610 -0,10943 0.02765 0,02072 -0,00638 

0,08155 -0,11651 0,02600 0,02052 -0.0049fl 

0,07465 -0,lJH4fi 0,01690 0,00604 -0,0174/ 

0,0/540 -0.12313 0,01074 0.01322 -0,02307 

0,07689 -0,1285b 0,02566 0,00754 -0,0144/ 

0,07513 -0, 12660 0,03202 0,01455 -0,0127/ 


6 

7 

8 

9 

10 

0.0953/ 

0,03633 

-0.00310 

-0,08512 

-0,03010 

0. 10088 

0,03138 

-0.06/31 

-0,07121 

-0,01960 

0,06914 

0,05441 

-0,00045 

-0,0/597 

-0,03118 

0,06098 

0,05933 

-0,09214 

-0,07871 

-0.03996 

0,05405 

0.05592 

-0,09427 

-0,06701 

-0,05194 

0,07786 

0.04830 

-0,11796 

-0,06664 

-0,05631 

0,06080 

0.03958 

-0,12597 

-0,04867 

-0,0338/ 

0,08584 

0,02693 

-0,09492 

-0,05020 

-0,04/09 

0,08869 

0,04096 

-0,07853 

-0,03351 

-0,04157 

0.07084 

0.01389 

-0,10210 

-0,04235 

-0,03442 

0,0fl7Q6 

0,00690 

-0,09406 

-0,06432 

-0,02599 

0,07332 

-0,00367 

-0,10047 

-0,06084 

-0,02637 

0,0824/ 

-0,00552 

-0,10223 

-0,07/50 

-0,02221 

0,08573 

-0,02953 

-0,09850 

-0,0784-3 

0,00520 

0,10622 

-0,02545 

-0.05576 

-0,08011 

0,01504 

0.10933 

-0,02603 

-0,04640 

-0,10906 

0,02331 

0,0830/ 

-0,02/14 

-0,04345 

-0,12389 

0,04189 

0,07260 

-0,05329 

-0.06697 

-0.12292 

0.06604 

0,08121 

-0,04615 

-0,04319 

-0,10829 

0,05742 

0,07923 

-0,03812 

-Q.05J75 

-0,09016 

0,0897/ 

0,07653 

-0,02952 

-0.03939 

-0.06727 

0,1 0382 

0,06653 

-0,02578 

-0,00529 

-0,03242 

0,09199 

0,03379 

-0,03870 

0,03123 

-0,00922 

0.09121 

0,00856 

-0,04879 

0,06195 

-0.00363 

0,07662 

-0.03800 

-0,02293 

0.08345 

0.01606 

0,07566 

-0,03499 

0,00461 

0,09586 

0,03369 

0,06826 

-0,03587 

0. 00523 

0,06766 

-0.01061 

0.04006 

-0,01573 

-0,01878 

0,06919 

-0,01382 

0,03762 

-0,03037 

-0,02066 

0,06908 

0.00347 

0,03729 

-0,04650 

-0,02317 

0,04850 

0,02323 

0,02160 

-0,06398 

-0,03100 

0,03560 

0,04521 

-0,01563 

-0,07835 

-0,03050 

0,02270 

0,07204 

-0,00043 

-0.08195 

-0,01168 

0,02576 

0,10139 

-0,01027 

-0,06011 

-0.04027 

0,03429 

0,14344 

-0,02347 

-0,0917b 

-0.03820 

0,03412 

0,13'304 

-0.02496 

-0.10606 

-0,05156 

0.050/4 

0.13612 

-0.01434 

-0,09475 

-0.04884 

0,0/472 

0,12276 

-O.Ol 142 

-0,06786 

-0.05169 

0,09184 

0,09102 

-0,06114 

«0,0o93l 

-0,01236 

0,1 0649 

0,10293 

-0,06703 

-0,08551 

-0.00506 

0,10177 

0,11537 

-0.06501 

-0,04904 

0.04181 

0,ll0Ot> 

0.07919 

-0.05/39 

-0,04790 

0,05229 

0,12020 

0.0/513 

-0,07473 

-0,04701 

0.06427 

0.10036 

0.09149 

-0,07256 

-0,071 40 

0.05738 

0.08451 

0,09260 

-0.07081 

-0,07080 

0.07563 

0.07306 

0.0919/ 

-0.04836 

-0,08988 

0,07529 

0.06} 17 

0,08185 

-0.022/3 

-0,08721 

0.08160 

0,06637 

0,07664 

-0.02960 

-0,07675 

0,06265 

0,0581/ 

0,06545 

0.00963 

-0,06889 

0.06956 

0.06113 

0,02/29 

0,00040 

-0,0540.5 

0.07677 

0,05149 

0,04233 

0,01602 
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TABLE 3 - LISTING OF THE FIRST 20 EIGENFUNCTIONS (continued) 


tlSkr^FUNCTJOKiS 11 - 20, POINTS SI • 100 I THfc V COMPONENT 


U 


la M 


EIGENFUNCTION 
IS u 


17 


IB 19 


0,0568ft 
0,06520 
0,06152 
0,06971 
0«05o6l 
0,02464 
0«U2664 
0,024^4 
0,02507 
0,01233 
• 0,00062 
0,00742 
0,03813 
0,04297 
Q,02492 
0,01223 
0,0257ft 
0,00131 
wO, 02767 
«<0,0i6S4 
-0,01423 
0,00fll9 
0,05514 
^O,O0S57 
!»0, 01126 
0,05626 
0,02295 
0,02250 
0,05366 
O,0ift27 
0,01151 
0 , 0020 ft 

-0,00676 

-0,04551 

-0,05922 

-0,03563 

-0,04570 

-0,02605 

-0,02404 

-0,06756 

-0,062S9 

#0,02755 

-0,01426 

0,01045 

-o,o1b69 

-0,01641 

-0,00931 

#0,02175 

-0,025B5 

-0,0224ft 


#0-03221 
-0,00759 
-0-0590ij 
0,00172 
0,00593 
-0,01005 
-0-02759 
*0,05509 
#0 , 0674fa 
-0, 09240 
-0,0ft«4B 
•0,05497 
#0,02B4b 
-0,01415 
0-00531 
C-0i26b 
0,02994 
0,064ft/ 
0,0794/ 
0,07203 
0 ,06461 
D, 08047 
0,03621 
0,00949 
0,03352 
0.0459fa 
0,0694j 
0.06393 
0 , 1 C 1 1 3 

0.08332 

0.059B2 

0,03039 

0,00212 

0,01BBb 

0,03/24 

0-01906 

#0,01747 

-0,04884 

#0,06615 

-0-08B9/ 

-0,08693 

-0,06202 

#0,09611 

-0,07735 

-u, 03219 
-0-0224U 
-0,00046 
0,02653 
0,0049/ 
0.03229 


-0.02521 

#0,00367 

0,02219 

0,01144 

0,02020 

0,00421 

0,00093 

0,00283 

0.03434 

0,03657 

0,03400 

0,03472 

0,01679 
0.02313 
O.OSOOfi 
0,01867 
0,07382 
0.065'15 
0,01470 
0,0 045*3 
-0,03287 
#0,06907 
-0,06153 
-0.06731 
-0,06087 
-0,0524ft 
-0,03801 
-0.03341 
#0,01362 
0,00864 
0,03047 
0,04355 
0,01003 
0,00444 
0,02829 
0,00221 
-0,00501 
-0,03122 
#0.01636 
#0.00717 
0,01669 
0.03972 
0,06217 
0,04)40 
0,01661 
0,01690 
0,00946 
0.00329 
0,00936 
0,00198 


-(5, 13979 
-0 ,14288 
-0,13783 
-0,13B24 
-0.11674 
-0,14603 
-0,10400 
#0,07151 
-0.C5409 
-0,05253 
-0,017/9 
-O.D) 12b 
0,01721 
0,03619 
0,05542 
0,06276 
0,07250 
0,09419 
0. 09660 
0.10054 
0,llB74 
0.00726 
0.09052 
0, 1 l59o 
0.14/59 
0,12490 
0,09385 
0,06535 
0,04915 
0,05/00 
0,04620 
0,06993 
0,06441 
0,01854 
-0,00562 
0,00340 
0,00345 
-0,04710 
-0,01530 
-0,02/54 
- 0,03260 
#0,02637 
0,00263 
-0^00866 
-0,05234 
-0,06993 
-0,12622 
-0,13094 
-0,1)317 
-0,0ftft73 


0, 0433ft 
0,01617 
0,01515 
0 , 00622 
0,02933 
0.00290 

#0,033B9 

-0.05447 

0,00065 

0.01107 

0,03291 

-0.00556 

-0,03466 
-0.03461 
0.02223 
0,02502 
-0.033H9 
-0,06405 
-0,05641 
-0.06107 
-0 ,04952 
-0.01 762 
0.07749 
0.12059 
0.13B56 
0-16224 
0.12634 
0,20796 
0,20636 
0,17941 
0.15614 
0.09279 
0,02961 
-0,02331 
-0,03874 
#0.04412 
#0.03601 
-0,07342 
#0.10108 
-0,12764 
-0,11 loO 
-0,08640 
#0,05643 
#0 ♦ 0 6fl85 
-0.07384 
-0,Qfi051 
#0,10361 
-0.06941 
# 0 , 0722 / 
-0,06243 


0,02371 
0,06357 
0,08155 
0,07683 
(U06fl0l 
0,04231 
0,0270 / 
0,02513 
0.00B09 
#0^00447 
-0, 01391 
#0,05653 

#0,()B90 4 
-0, 08016 
-0,05483 
#0.08287 
-0,1348/ 
#0.15062 
#0,13311 
-0,15579 
-0,13729 
-0.06570 
-0,03601 
0, 0201ft 

0.06844 
G, UO36 
0,1091/ 
0, U94i 
0,11039 
0.05303 
0,07655 
0,11746 
0.12764 
0,11815 

0,07579 

0,00393 

-0,0584/ 

#0,07669 

#0,01999 

-0,05176 

-0,05443 

#0,02204 

0,01263 

#0,02326 

-0,01734 

0,00653 

-0,00596 

-0.00143 

-0,01829 

#0,02042 


-0.03005^ 

-0,02942 
0.0242} 
0, 03872 
0,02446 
#0.04136 
#0,04200 
#0,01748 
0,00643 
-0,04376 
#0,02906 
-0,02642 
-0.01520 
0.00141 
0,02456 
0,04990 
0,00869 
-0.Q2412 
-0.02/42 
-0.01386 
0,01248 
-0,00829 
#0,0J45B 
-0,01465 
0.04730 
0,006B9 
#0,03659 
-0,03310 
-0.00470 
0,04247 
-0.00597 
0,0l90fl 
0.05297 
0,05681 
0,06945 
0.06567 
0,07190 
0,00072 
0,01623 
#0,0 0355 
-0,00139 
0,03483 
0,01147 
-0,03579 
#0,07152 
-0.05015 
0,00091 
-0,01920 
#0,055/1 
-0.05467 


0,00472 
0,02092 
0^09733 
0, 09456 
0,09594 
0,07557 
0.05622, 
0,01647 
0,03727 
0,04586 
#0.01011 
-0,04293 
-0,03304 
-0.03179 
-0,03136 
#0,05/67 
-0,04250 
-0,06223 
-0.0534B 
-0,04761 
-0,01705 
0.00056 
0,021SB 
0,0ft216 
0,12635 
0,13655 
0, ].539ft 
0, 08314 
0,06067 
0,02191 
0,02900 
-0,01B67 
-0 , 05556 
-0.09414 
-0,09609 
#0,06779 
#0,07033 
-0,04676 
-0,03461 
#0,06596 
-0,06284 
-0,01107 
-0,00101 
-0,01 3B2 
-0,04/97 
#0.06537 
#0,05120 
-0.04695 
-0.05407 
-0,02B16 


0,04154 
0,0135B 
0,02345 
-0,01912 
#0,03452 
#0,03602 
0,01)041 
-0. 01940 
0,01396 
0,02862 
0,03750 
0 ,00559 
0,01011 
0,00603 
-0.QQ695 
-0,01063 
-0.04375 
-(1,03486 
-0,01326 
-0, 02691 
0.00263 
-0,02451 
#0,07442 
-0,07348 
#0,09522 
#0 ,05606 
0,00616 
0.05253 
0.01199 
-0,04229 
-0,0/679 
-0,06888 
-0,12476 
-0,090/0 
-0,06560 
-0,0/556 
-0,04132 
-0,04050 
(),00 J 36 
0,04739 
0,00/44 
0,000/5 
0,08463 
0,06599 
0,06792 
0,09109 
0,10069 
0, 10969 
0.06/ 03 
0.06832 


0 

-0 

#0 

-0 

-0 

#0 

-0 

-0 

0 

0 

0 

0 

0 

0 

a 

0 

0 

0 

0 

-0 
— 0 

0 

0 

0 

0 

0 
-0 
-0 
-0, 
-0. 
#0. 
-0, 
-0. 

# 0 a 

# 0 , 

-Oa 

o< 

0 . 

0 . 

01 
Oi 
0 . 
0 . 
Oa 

- 0 . 
# 0 , 
- 0 . 
-Oa 
— 0 a 
-Oa 


20 

,01141 

,04162 

,03354 

,00746 

,05353 

,05983 

, 03/22 

,00382 

,03761 

,0}339 

,02224 

.Qlbtl 

.00913 

,03891 

,04002 

,02305 

,02719 

.00529 

>00046 

,0422/ 

,05702 

,00009 

^3151 

,0764/ 

,09524 

.04841 

,07603 

,06249 

>04502 

.04H61 

,06350 

,0405ft 

,03539 

,04572 

,Q3l2i 

, 01/22 

,01649 

.08741 

,09193 

,06622 

.02757 

,06258 

,0615ft 

.03471 

,00367 

,02365 

,02461 

.03324 

,04011 

,06932 





1 

0t0B600 

OpOB'Jab 
0 ,0BH2b 
0,08‘?02 
0,0^047 
0.0^276 
0,09660 
0. 100B9 
0.10^08 
0,10679 
0 . 106^1 
0,10954 
0,11524 
0,10907 
0,10764 
0,10999 
0,10875 
0,11215 
0, lOBSi 
0.11732 
0,11530 
0,10745 
0,10546 
0,09849 
0,08851 
0,09134 
0,10443 
0,10760 
0,11196 
0.11382 
0,11688 
0, 10V34 
0,10858 
0,10710 
0.11014 
0,11510 
0.12037 
0,11868 
0.11854 
0,12790 
0.12178 
0, llBBg 
0,10921 
0.10326 
0.09999 
0.09878 
0,10490 
0, lOlSH 
0,10263 
0,09859 


TABLE 3 “ LISTING OF THE FIRST 20 EIGENFUNCTIONS (continued) 


LISENFUNCTIQNS I - 10, POINTS 1 - 50 I TH£ U COHPU^fcNT 




EJBENFUNCTION 





2 

3 4 

5 

6 

7 

6 

9 

10 

-0.07741 

0,07562 0.0B443 

<-0,05637 

"-0, 09798 

<*0,05196 

-0,10021 

-0.08565 

-0,06500 

-0,07105 

O.OB128 0,0769b 

-0.05677 

-0.09849 

-0.08333 

-0,09487 

-0,08679 

-0,07833 

-0,07269 

0,08379 0,06700 

-0.07093 

-0,09637 

-0,08531 

-0.04756 

-0,11991 

-0,06447 

-0,08000 

0.06512 0,07255 

-0.06657 

-0, 09813 

-0,07565 

-0, U3345 

-Q.U229 

-0 • 06650 

-0,06799 

0.07999 0,0773b 

-0.07871 

^0,10165 

-0.07576 

-0,03736 

-0,10356 

-0,08196 

-0.08312 

0.07994 0.08178 

-0,07279 

-0,09125 

^0,06434 

-0,06110 

-0,08650 

-0,07537 

-0,08054 

0,08384 0,08721 

-0.06469 

wO, 0/929 

•<0,05539 

-0,08143 

-0,08668 

-0.Q4603 

-0, 0640b 

0,08001 0.0808B 

-0.07941 

-0,06BB6 

-0,05295 

-0,10641 

-0,04376 

-0,05933 

-0,08485 

0,07719 0,07085 

-0,09543 

<•0,0 4663 

-0,07084 

-0,09563 

(},00252 

-0,08265 

-0.08612 

0,06979 0,07173 

-0.P9797 

-0,02975 

-0,05719 

-0,11337 

0,02540 

-0,06612 

-0,09369 

0,07081 0,07404 

-0.09570 

^0,04779 

-0,06157 

-0,11630 

0,03681 

«0«06224 

-0.09459 

0,06154 0,079/7 

- 0,07743 

-0,05249 

-0,06316 

-0,09220 

0,02424 

-0,07262 

-0,08652 

0,07078 0,06421 

-0.07402 

-0,04565 

-0,06416 

-0,10036 

0,02464 

-0.05456 

-0,08964 

0,07060 0,06163 

-0.08554 

-0,02565 

-0 .04665 

-0.10742 

0,02933 

-0,03600 

-0.08916 

0,06693 0.04358 

-0,07158 

^0,00854 

-0,01785 

-0,11010 

0,02685 

0,00167 

-0,09124 

0,06020 0,02572 

-0,07183 

0,0169/ 

-0.01553 

-0,08331 

0,03912 

0,03343 

-0.08/31 

0,05683 0,03574 

-0.09661 

0,02472 

-0,00/02 

-0,05105 

0,03856 

0,05160 

-0,08866 

0,06516 0,03873 

-0,08496 

0,03394 

0.01619 

-0,06197 

0,02636 

0,10008 

-0,08656 

0,06043 0,04758 

-0.07370 

0,04299 

0.04519 

-0.02786 

0,02300 

0,12699 

-0.06'/4i 

0,05127 0,03416 

-0.07072 

0,02443 

0,04313 

0,00529 

0,01686 

0.13764 

-0.0770i 

0,04539 0,02757 

-0,08185 

0,00410 

0,04667 

0.03075 

0,01889 

0.15187 

-0.07519 

0,03119 0,02256 

-0,08776 

0,02672 

0,06727 

0,01022 

0,03455 

0,19174 

-0,07225 

0,03647 0,02665 

-0,08153 

0,03543 

0,06796 

0,01576 

0,03416 

0,16943 

-0.0631 6 

0,02665 0,01585 

- 0,04939 

0,0751 1 

0,05051 

0,01519 

0,06945 

0,15511 

-0,05040 

0,01134 0,00862 

-0,05730 

0,08921 

0,06897 

0,02B53 

0,o53) 3 

0,12156 

-0,04964 

0,02568 -0,03228 

-0,06250 

0,09066 

0,05052 

0,02550 

0,01480 

0,09192 

-0,05464 

0.01553 -0,04167 

— 0 . 02866 

0,10413 

0.01914 

0,01427 

-0,018U7 

0,08504 

-0,05985 

0,03069 -0,04648 

-0.00043 

0,09601 

0,02718 

0,01798 

0,02455 

0,13679 

-0,06945 

0,02059 -0,06005 

0.01568 

0.1 01 73 

0,04007 

0,02505 

0,0397/ 

0,16896 

-0,0636b 

0.01767 -0,06072 

-0.00836 

0,07841 

0.03108 

0.03368 

0,0608b 

0,12366 

-0,06333 

0,00234 -0,07567 

-0.00878 

0,045RB 

0.03471 

0,05550 

0.03848 

0,10503 

-0,0621 1 

0.OD2B3 -0,06241 

-0,01611 

0,051 /9 

0,04522 

0,09601 

0.05634 

0.10989 

-0,06273 

0,01188 -0,07934 

-0, 01815 

0,01251 

0,0369a 

0,09472 

0,04859 

0.08735 

-0, 06323 

0,02403 -0,09309 

-0,02229 

-0,00043 

0,04322 

0,09733 

0,03296 

0,06369 

-0,06106 

0,03699 -0,08913 

0.00832 

0.00775 

0,03985 

0,10242 

0,03948 

0,01201 

-0.05630 

0,04363 - 0,10643 

0,02739 

0,01933 

0,06785 

0,11506 

-0,00341 

0,00042 

-0,05531 

0,05674 -0.104B1 

0,03664 

0,02868 

0,04759 

0,11417 

-0,01051 

0.00596 

-0,04659 

0,05354 -0,10408 

0.03295 

0,03704 

0,04548 

0,09121 

-0,02262 

0,00044 

-0.04407 

0,05399 -C, n 105 

0.02865 

0,00235 

0.04747 

0,07894 

-0,04862 

-0,04187 

-0,0384(; 

0,05439 -0,12490 

0,03538 

-0,01061 

0,05389 

0,0547? 

-0,08416 

-0.06058 

“0.0413 fa 

0,04969 -0,13649 

0.06071 

-0,01653 

0,04122 

0,05691 

-0,06674 

-0.10723 

-0.04272 

0,05244 -0,14385 

0.06101 

-0,02349 

0,02634 

0,06036 

-0,01416 

-Q. 06614 

-0.04435 

0,05976 -0,16532 

0,06047 

-O.0199U 

0,02594 

0, 05867 

-0,01742 

-0,11622 

-0.05186 

iJ, 06349 -0,1 5429 

0. 0R42O 

-0,01 829 

0,02744 

0,05787 

-0,02612 

-0,13093 

-0.05373 

0,06435 -0,14023 

0.09616 

-0,01250 

0.01218 

0,02777 

-0,01986 

-0.1 2868 

-0.06171 

0.05467 -0,14844 

0.11655 

0, 01812 

0,00051 

0,02991 

-0,00796 

-0.14333 

-0,05489 

0, 1)3309 - 0.13795 

0.12401 

0,01224 

-0,01774 

0,03802 

0,00073 

-0.15082 

-0.05/66 

(>,03882 -0.15488 

0. 1 1683 

0,00752 

-0,02258 

-0,00517 

-0.00251 

-0.12184 

-0.05B26 

0,03628 -0,14094 

0.09593 

0,0U83 

-0,02699 

0,00433 

0,00676 

-0.12253 

-0,05/43 

0.04362 -0,13367 

0.07950 

0.0132/ 

-0,01306 

0,01193 

0,01513 

-0,14284 
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TABLE 3 - LISTING OF THE FIRST 20 EIGENFUNCTIONS (concluded) 

tlufcNFUNCTIOr-'S H - 20 , POINTS lOl - 150 | THE W COMPONENT 

ElfiENFUNCTlOfl 


n 12 13 14 15 

0,10667 -0.19B81 0.1«U9 0,18613 -0.08965 

-0,04926 -0.16964 0,14633 0.124B3 0,0819/ 

-0,04917 -0,05649 0,12257 0,03/69 0.06482 

-0,02942 -0.03542 0,014/(6 -0,fr3O29 0.03893 

•0,0(?980 0.C3O/4 -0,061/10 -0,04942 0,08115 

-0.023/9 0.12501 -0,08023 -0,0/718 0.06085 

-0,03133 0.17024 -0,05585 -0.05120 0,03160 

-0,01752 0.18222 -0,08415 -0,04322 -0.03559 

-0,08537 0.U8068 -0.13553 0,00l30 -0,10/52 

-0,04823 0,09770 -0.21427 0,02445 -0,19086 

-0,02950 -0,02743 -Q. 19884 -0.016/7 -0.10962 
-0,01618 -0.16612 -0,04662 -0,0/068 -0,02464 
-0.03529 -0.143/8 -0.05334 -0.04039 -0.04160 
-0,02803 -0,193/2 -0,05500 -0,04024 0.023B6 

0, 01049 -0,17560 -0,05.587 -0.07295 0,0/985 

0,08268 -0,11797 -0.02150 -0.0j2/i 0,16152 

0, 03163 -0,1479a 0,03026 -0,00727 0.08499 

0,08644 -0,07232 0,01766 -0.01042 -0.02401 

-0,00171 0,06301 0,03376 -O.0B789 -0,10070 

-0,0l38l 0,l307o 0,14226 -0.J2695 -0,06356 

0,0/304 Q, 08089 0,14571 -fi. 01969 0.006/5 

0,07570 0,12262 0,09949 0,038/3 «0,025}6 

0,13536 0,15553 0,04180 -0,12669 -0,03958 

0,11576 0,15821 0.00046 -0,14055 -0,00271 

-0.00598 -0,08171 -0,07962 0,06572 -0.00695 

-0.01026 -O.0839i -0,11978 0.2026/ -0.03400 

0,01222 -0.10130 -0,00064 0,06219 0,06845 

0,01869 „0, 08366 0,10464 -0,14058 0,15882 

0,03224 0,09013 0,16631 -0.20828 0.23190 

-0,02359 0,15600 0,14696 -0,20445 0.12139 

-0,16921 0,231 03 0.105('3 -0,04664 0,1 2175 

-0,28446 0.13273 0,17B()0 -0,06099 o.065o/ 

-0,2BU9 0,02249 0,j 6905 0,04/28 -0.09589 

-0,23209 -0.149ms 0,139i) 9 0,in2t^2 0.05469 

-0,19207 -O.ibbbo 0,10041 0.17593 0.07153 

-0. 19272 -0.04224 0,04859 0.11923 -0.02541 

-0,04584 0,03539 -0,12566 0,l0559 -0,07590 

0, 05734 0,02209 -0.15821 0,03299 -0.10953 

0,02390 -0,08731 -0,19068 -0,07414 -0,l0344 
0,01105 -0,05687 -0,19780 -O.U700 -0.10231 
0,03695 -0,05063 -0.27091 -0,05206 0,02192 

-0,01756 -0,03135 -0,21.497 -0,04779 0.07466 

0, 02313 0.01880 -0,08095 -0,00028 O.UflOO 

0,11653 -0,02387 -0,01951 0.02666 0,09252 

0,20281 -0.05299 0,02618 0, 0666b 0,05107 

0,22151 -0,06731 0,09186 0.06/34 0,00690 

0,22464 -0,03583 0,12392 0,05893 -0.09.449 

0,16157 -0^07/57 0,14140 0.00105 -0.19731 

0,193.36 -0.10422 .0,17120 0,00997 -0.1741 1 

0,1/056 -Q,10204 0,18126 0,03662 -0.05179 

0,01031 -fO, 05298 -0,06590 -0,04472 0,05954 


16 17 18 19 20 

-0,09415 0,01297 0.02248 -0,10424 -0,07110 

-0,07028 0,16882 0,U868 -0,18638 0,01249 

0,01106 0,11969 0,06269 -0,12641 0,0960/ 

-0,06073 0,00575 0,00725 -0,11335 0,1861/ 

-0,06238 -0,03922 srO, 14716 0.02531 mi),QQUS 

0,05133 ^0,08046 -o;o4037 0,09259 0,OQ997 

0,0740ft -0,13642 0.03901 0,0/030 0,11485 

0,12860 -0,15715 -0,06108 0,01501 0,16564 

0,08541 -0,05750 -0,0/720 -0,00260 -0,01270 
0,02/94 0,06086 •^0.03515 -0,12467 -0,15774 

-?0, 01456 0,11363 0,00432 -0,03607 -0,12646 

0,03496 0.22602 -0,04655 0,01090 -sO. 12945 

-0.04526 0,26636 0,01912 0,08620 -0,07486 

•^0,011 90 0,14719 0,19345 0,l3Sil 0,08201 

-0,03065 -0.03388 0.10232 0,14555 0,10128 

-0,02554 -0,11088-0,01720 0,25246 0.12797 

0,01/21 ?-0,2l7SB -0,02976 0,23240 -0,05334 

-0,02383 -0,27231 0,06290 0,06287 -0, 11522 

-0,04368 -0,2283? 0,14493 -0,2/349 -0,l09l5 

■!*0. 02672 -0,10762 0,16261 -0,1/154 0,08/09 

-0,12117 -0.06906 0,10525 -0.11578 0,11476 

-0,0755X 0,13063 -0,11046 0,04790 ..0,02055 

0,18546 0,27644 -0.14905 0,06026 -O.OlOli 

0,10364 0,19544 -0.15590 0,04249 -0,15444 

-0,01889 -0,01852 -0,04/89 0,08797 -0,04621 

O.OSOpi -0,00099 0.05151 0,06137 0,09699 

0,07.301 -0,1 1 01 4 0,01952 -0,19862 0,10997 

-0,16954 -0,21684 -0,0/046 -0,21904 -0,07494 
-0.18464 -0,05281 -0,01514 -0,12775 -0,19740 
-0,05/82 0,168/4 -0,03b52 -0,05243 -0.01226 

-0,0010/ 0,19001 0.11122 0,12898 0,18935 

-n,05/0'J 0.09103 0,08564 0,13503 0^08918 

0,00838 -0,04184 0,02130 0,05325 0,0638/ 

0,079o5 -0.17952 -U.l39o0 0,10/96 0,08080 

0,11443 -0.12190 -0.14807 -0,05114 -0,09252 
0,06301 -0,00395 -0,2293/ -0,04219 -0.06196 
0,13538 -0.01949 -0,12290 -0,13601 0,14224 

0.18551 -0.0791b 0,0351 1 0,00278 0,13.388 

0,05543 0,00678 0,18611 -0,05560 0,00677 

-0,14165 -0,02905 0.29623 -0,12250 0,03'/5l 

-0,20601 0.08363 0,06454 0,05082 -0,03863 

-0,30438 0,06994 -0,01712 0,03726 -0,02406 

-0,19666 -0,06356 -U, 11390 0,03212 -0,12499 

-0,0356/ -0,13460 -0,00756 0,06021 -0,26212 

0,12880 -0.058/2 0,00407 0,11633 -0,17709 

0,02738 0,00664 0,15550 0,12561 -0,13030 

0,00002 0,11362 0,14139 0,00255 0,0l94/ 

0,01480 0,10296 0,00281 -0,09223 0,10125 

n,103}5 0.04366 -0,00560 -0,02390 0,02443 

0,1097/ 0,0i454 -0,09021 -0.06997 -0.04999 

-0,10501 -0,03323 -0,14053 0,02210 O.OOOOO 
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5.2 The Sampling Properties of the Expansion Coefficients 

Forming a 100 150 matrix F of the 100 ensemble functions and a 

150 20 matrix E of the 20 first eigenfunctions and performing the 

matrix multiplication 

F X E = B 

gives us the 100 x 20 matrix, B, consisting of the expansion coefficients, 

a^. A histogram is then calculated for each of the 20 columns and the 

corresponding probability density function estimated. The goodness of 

the fit can be tested by various tests. 

In this analysis, a subjective study of the histograms revealed 

Gaussian probability density functions for the expansion coefficients, 

so that the only statistics necessary are the means and the variance 

of the 20 columns. If we had prepared the ensemble in such a way that 

E{f } = 0, we would have had 
n 


E{a } = 0 
n 


(5.2) 


E{oc a, } = X 6 - 
n k n nk 


and due to the Gaussian assumption all the probability density functions 

would be known, in a statistical sense, because A , which equals the 

n 

til 

variance of the n expansion coefficients, is known. In this analysis, 
the relation (5.2) provides an internal check on the calculations because 
for one expansion coefficient, the variance plus the squared mean must 
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eqiial the corresponding eigenvalue. Table 2 gives the mean and the 
variance for the 20 expansion coefficients. As would be expected, 
only the few first have a mean significantly different from zero. 

The mean and the variance for a particular expansion coefficient 
are only estimates of the true values, and by standard statistical 
methods we can obtain 95% confidence limits for both quantities; in the 
application of the model we can investigate the effects of choosing 
values inside as well as outside these limits. 

The apparent Gaussian distribution of the expansion coefficients 
needs further discussion, especially since Smith (1971) reports he found 
them to be uniformly distributed. 

When the domain of definition is extended to and R(s,t) *= 

R(s-t) , then (4.3) becomes 

00 

[ R(s~t) (J) (s) ds = A c() (t) 

J ^n n^n 

— OO 

but in this limit the eigenvalue spectrum becomes continuous, and we 
can write 

00 

I R(s~t) (j)((A),t) dt = A(w) 4>(a), s) 

--00 

with the solution 4)(oj, t) = e Hence, A(aj) is the Fourier transform 

of the correlation function, the usual power spectrum. And so the 
expansion coefficients become in the limit equal to the usual Fourier 
coefficients, which for all practical purposes can be assumed to have a 


Gaussian distribution. 
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5.3 The Probability of Occurrence 

A direct assessment of the probability of occurence of generated 
time sequences which are critical to the system being studied is not a 
straightforward procedure with the model developed here. 

However, in order to analyze the problems let us consider the 
theoretically simpler model in which an ensemble of k records each of 
N points, k ^ N, with the correlation function R(s,t) has led to N 
eigenfunctions each specified at N points and to N probability 
density functions, one for each of the N expansion coefficients 
°^1* **** ■ Artificial time histories are generated using the 

eigenfunctions and by sampling in the probability density functions. 

The generated series will also have the correlation function R(s,t), 
and the first probability to be determined is the probability that the 
systems we are concerned with will encounter turbulence with the given 
R(s,t). This probability can presumably be estimated by considering 
the measured ensemble in relation to the total collection of empirical 
data about turbulence. 

The next probability we must determine is the probability of ob- 
taining a critical time sequence by sampling in an ensemble of functions 
with the correlation function R(s,t). 

Let us assume that all the critical functions with this correlation 
function are in the subspace of the total phase space where 

1 “l i ajj < Ojj < + Aajj) 

(-00 < < 00 , -00 < < oo ) 
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Here Prob(a^,a 2 > ... e fl) = 1, and the probability we want to 

determine is Prob (ot_ , a«, ... e In order to limit the 

X z JN A 

search for we realize that we are only interested in probabilities 

that exceed a certain limit, P„ say. The part of in which the 

t* 

higher order expansion coefficients differ significantly from zero 
is very likely to be associated with very small probabilities (see 
Table 2) , hence it is justifiable to search for only in the space 

% “2* “m* ^m+1 ^ ...» = 0) 


where the limits for a_ , ...» ot are connected with P„ . 

1* m E 

Suppose that in addition to the lack of correlation expressed by 
(5.2), the coefficients are also Gaussian distributed, an assumption 
justified in Section 5.2. Then the coefficients are independent and 
after having established by trial and error the boundaries of the 
probability we are searching can then be estimated from the Gaussian 

P, 


probability density functions p^, p^. 


Prob (a-, a«, a £ fi.) = 

12 m A 


m 


a +Aa 
m m 


Pi d«i 


p2 dOi . . . 


p da 
m m 


m 


Mainly because of computer limitations, we are forced to generate the 
time-sequences in small pieces which are patched together and then to 
perform transformations in Fourier space in order to move energy to low 
frequencies. Hence it is obvious that the above described procedure for 
obtaining the probability of occurence of critical time sequences cannot 
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be applied to the model formed In the preceding sections. We will 
have to rely on physical arguments, not as much to calculate proba- 
bilities, but merely to establish whether the generated critical time 
series appears to be physically realizable. 
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6.0 DESCRIPTION OF THE TURBULENCE RECORD 


The turbulence record was obtained by analog instrumentation during 
the 1969 Kansas Experiment by N. E. Busch and S. E. Larsen (Larsen (1971), 
Busch and Larsen (1972)). It was measured at 5.66 meters height 
(between 15:31 - 16:31 on 30 July, 1968), The stratification was 
slightly unstable, the Richardson number being -0.101. The mean wind 

speed was 6.56 m/sec and the variances for the fluctuations u, v, w, and 

2 -2 2 -2 2 “2 
T (temperature) were 2.414 m sec j 1.904 m sec ; 0.359 m sec- ; 

0.761® C^. 

The analogue signal was later digitized at 1000 Hz and transferred 
to digital tape. For the purpose of this analysis, the signal was 
further block averaged over 100 points to give a 10 Hz signal. A 
sequence of 50 minutes was selected giving 30,000 data points for each 
of the components u, v, and w (the temperature data was not used in this 
study). Figure 6.7 shows the whole record with 1 second block 
averaged values, and Figure 6,8 shows 100 seconds of the record with 
10 Hz values. In the last figure, all the data has been normalized 
to mean zero and variance 1. The analogue signal was reversed In time 
during the handling and so we will expect our model to generate time re- 
versed turbulence. A simulated record should then be time reversed 
before applied in practice. 

In order to assess whether the generated turbulence behaves like 
real turbulence or not, a set of criteria were proposed in section 1. 

One of the criteria requires the model to produce signals that possess 
the notable observed statistical characteristics of observed turbulence. 
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In the rest of this section some of these chai 'cteristics will be 
estimated from the turbulence record and the t ■ . f5 will be compared 
with those obtained from a similar analysis by Dutton and Deaven (1971) 
(hereafter referred to as DD) . 

6.1 Probability Densities and Distributions 

The use of the probability functions is extended to powers of the 
velocity variables, here to the fourth order. Figures 6.1 and 6.2 
show the probability density and distribution functions for the 
standardized velocity components. In each graph, the solid line 
illustrates the Gaussian behavior with the curves for higher orders 
derived from the transformation 



DD reached the conclusion that the frequency functions departed 
from those of a Gaussian process from a set of figures which showed a 
behavior quite similar to Figures 6.1 and 6.2. 

6.2 Increments 

Using one of the component series, for example u, we can construct 
a new series 


( 6 . 1 ) 


u(x,L) = u(x) - u(x + L) 



PROBABILITY DENSITY 





Figure 6.1 Kansas turbulence. Probability density and 

distribution functions for the first and second 
power of the standardized data. The Gaussian 
case is illustrated by a solid line. Integers oh 
logarithmic vertical axis denote power of 10 



PROBABILITY DENSITY PROBABILITY DENSITY 
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D 6^ 12? n2 296 Q 6^ 12? 1S2 296 
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Figure 6.2 Kansas turbulence. Probability density and 

distribution functions for the third and fourth 
power of the standardized data. The Gaussian 
case is illustrated by a solid line. The integers 
on the vertical axis denote power of 10 


r 
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by means of which we can analyze spatial variations in the turbulence. 

The timestep At = 0.1 sec is transformed into a length step Ax using 
Taylor’s hypothesis: Ax “ uAt = 6.56 m sec ^0.1 sec = 0,656 meters. 

The structure function. Figure 6.6, is defined with these 
increments as 

(6.2) D(L) = E{u^(x,L)} 

and it should be dependent only on internal conditions of the flow. 

The statistical properties of the increments are Important in 
several ways. First they characterize the spatial variations in the 
turbulence- -hence the distribution of increments provides insight into 
the uniformity of spatial structure of the velocity fields; second, if 
the series u(x) are Gaussian, the distribution of Increments would also 
be Gaussian; and third, they are useful in the study of whether the 
data is self-similar in the sense of Mandelbrot as described in 
Section 3. 

The probability functions for the increments are shown in Figure 6.3 
and Figure 6.4 for lags 0.66, 6.6, 66, and 666 meters. Again, the curves 
are very similar to those given by DD, demonstrating that the increments 
with small lags deviate more sharply from Gaussian behavior than the 
Increments at large lags. 

6.3 Measures of Intermittency 

As argued by DD, the non-Gaussian behavior appears to be intimately 
related with the intermittency of the turbulence and they provided fur- 
ther measures of the intermittency. Dutton, Lane, et al. (1969) define 
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INCREMENTS L- ID M INCREMENTS L= ID M 


Figure 6.3 Kansas turbulence. Probability density and 
distribution functions for the standardized 
data. The intervals 1 and 10 m indicated on 
the graphs signify actual lags of 0,66 and 
6.6 m. The Gaussian case is illustrated by a 
solid line 
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an Intemiittent process as a process where a relatively large fraction 

of the variance is contributed by a relatively small fraction of the 

total record. DD then considered how moments such as the variance, 

skewness, and kurtosis accumulate as a function of the fraction of 

the total record length. These statistics are shown in Figure 6.5. The 

numerical algorithm can be thought of as a process in which the record 

is rearranged so that the observations are ordered by size; the 

curves are then obtained by summing the appropriate power of these 

observations and plotting the result against the fraction of the 

observations used in the sum, DD showed that the same curves can be 

obtained from the probability density function as follows: 

til 

The fractioii of the 2m moment (m = l,2,...n) contributed by 
observations with absolute value greater than jy | is 


- y 


-}-oo 




Py(x)dx + 


^2m p^(x)dx] 


^2m^y>= 


|yl 


r| 2m , . ,-l 

[ X Py(x)dx] 


and the fraction of record occupied by observations with absolute value 
greater than |y] is 


-|y 






p (x)dx + ' 

y 


yl 


p (x)dx 

y 


For odd .moments we have the contribution by observations with 
values less than y 
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y CO 

” J Py(x)dx [ I Py(x)dx] 

—CO —CO 

—CO 

Thus, if the density functions Py(x) are known, then both F and 
R can be determined as function of y and so F is known as function 
of R. 

Also, these figures are very similar to those given by DD with 
only a slight departure from the Gaussian case for the variance and 
a distinct deviation for the skewness and the kurtosls. The difference 
found between the components can also be found in some of DD*s figures. 

6.4 Exceedance Statistics 

Among the various exceedance statistics that can be used for 
studying statistical structures, DD choose N(y)/N(0) which is the 
ratio of the number of crossings of value y with positive slope per 
unit time to the number of crossings of zero with positive slope. 

Figure 6.5 shows N(y)/N(0), and again we find the same behavior as 
found by DD. 

6.5 Spectra and Covariance Functions 

The analysis of turbulence relies heavily on the theory of second 
order processes, the covariance function being one of the most inqtortant 
characteristics of such processes. From the covariance function, a 
generalized spectral representation can be obtained, as shown in 
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Figure 6.5 Kansas turbulence. Accumulation of variance, 
skewness, and kurtosls (on the vertical axis) 
as a function of the percentage of the record. 
Note that the final S 3 nobol on the right gives 
the average variance, skewness, and kurtosls 
respectively. Lower right graph: the ratio 

for the standardized data of the number, N(y), 
of crossing of level y with positive slope to 
the number, N(0), of crossing of zero with 
positive slope as a function of y/o 
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Section 4, where the shape of the eigenfunctions can be argued to 

be of some significance. Further, the eigenfunctions gave a unique 

representation of the correlation function, and in the limit where we 

could assume stationarity and infinite integration limits, this 

representation produced the well-known fact that the correlation function 

and the spectrum constitute a Fourier transform pair. 

By means of correlation fimctions and spectra, the analyst is able 

to investigate the sequence of events throughout an enormous amount of 

data by looking at smooth curves. Although the correlation functions 

and the spectra contain the same amount of information, the spectra 

are normally best suited for a subjective analysis because they reveal 

how the variance (or the covariance) is distributed over wave numbers 

and hence how the energy is distributed on scales. 

The autocorrelation functions R , R , R and the cross- 

uu w ww 

correlation functions R , R , R are shown in Figure 6.6 on a 

uv uw vw 

logaritmic lag scale. As found in DD, the correlation of the horizontal 
components u and v is higher in the midrange than that of the vertical 
component w. The R^^ function shows the expected behavior to tend 
to a negative value significantly different from zero at small lags 
indicating a downward transport of horizontal momentum (the Reynolds 
stress). Also shown on the figure is the structure function D(r), 
where for a spatially homogeneous process, the relation between R(r) 
and D(r) is easily found to be 

(6.3) D(r) = 2a^ (1 - R(r)) 

[ 

a relationship that can be seen to hold well for the curves in Figure 6.6 


L 
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Figure 6.6 Kansas turbulence. Autocorrelation, cross- 
correlation and structure functions for the 
standardized data plotted as a function of the 
lag on a logarithmic scale. The plotted 1 
refers, to the u-v correlation functions, the 
2 to u-w, the 3 to v-w. The integers on the 
horizontal axis denote powers of 10 
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The spectra S . S , S and the cospectra S are shown in 
u V w uw 

Figure 6.7. The most pronounced characteristic of the spectra is the 
--5/3 slope exhibited by the u and v components over one decade of 
frequencies. The w component is seen to flatten out at low frequencies. 
The cospectrum gives the wavenmnber decomposition of the Reynolds 
stress responsible for the transformation of mechanical energy and hence 
is an important statistic to model correctly.. The cospectrum is 
seen to be significantly different from zero at frequencies 0.01 - 1, 
indicating the active scales in the downward transport of momentum 
to be of the order 600 m to 6 m. 
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Figure 6.7a Kansas turbulence. First third of record; 

not normalized. One second averages plotted 
versus time; the numbers on the horizontal 
and vertical axis denote seconds and meters 
per second respectively 
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7.0 TESTING AND FURTHER DEVELOPMENT OF THE MODEL 


Two experiments with the model are discussed in this section. The 
generation of the turbulence follows the generation scheme given in 
Section 5.3. 

As a first step we choose to generate the series in sequences 
of length 30 sec (= 300 datapoints). This is six times the length of 
one of the component parts of the eigenfunctions and our choice is 
obviously motivated by the way the ensemble was constructed. In^- 

side every 30 second interval we construct the turbulence in pieces of 
5 seconds. One of the pieces, selected at random, is created as active 
turbulence and the five other pieces as passive turbulence. 

The differences between the two experiments are as follows : 

1. The uniformly distributed random series consisting of integers 
between 1 and 6 which give the position of the 5 sec active turbulence 
piece inside the 30 sec sequence is sampled for each experiment. 

2. A Gaussian distributed random series is sampled to obtain 
the expansion coefficients for each experiment. 

3. In experiment 1, the mean and the standard deviation for the 
expansion coefficient distributions are the same as given in Table 1. 

In experiment 2, the standard deviations were changed to '^2* ^3^ ~ 

(8.0, 8.0, 4.0). 

4. The alteration of the expansion coefficients in the five 
passive intervals was done by multiplying each of the three first 
expansion coefficients by a uniformly distributed random number 
between 0 and 1 for experiment 1 and 0.5 and 1 for experiment 2. 
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In both experiments, the third expansion coefficient was, for 
every 5 sec generated interval, multiplied by +1 or *^1 picked at 
random. The justification for this procedure lies in the shape of the 
eigenfunctions and the distributions of the expansion coefficients. It 
follows from Table 2 that in contrast to all the other expansion 
coefficients has a much larger mean than standard deviation, implying 
much greater probabilities for obtaining positive rather than negative 
values. Figure 5.1 shows that most of the contribution of the third 
eigenfunction is the peak in the w-component, which in creating of the 
ensemble was chosen positive. But we want to generate negative 

gusts as well as positive, and therefore the + 1 multiplication of 
was introduced. 

7.1 Statistical and Sequential Characteristics of the Two Experiments 
Some 8200 datapoints were generated for each component for each 
experiment and then subjected to the analysis applied to the turbulent 
record in Section 6. The result is displayed in the Figures 7.1 ~ 7.18. 

The probability densities and distributions in Figures 7.1 - 7.14 
show no large deviations from experiment 1 to experiment 2 and a 
comparison with the Figures 6.1 and 6.2 reveals the generated turbulence 
to exhibit the expected non-Gausslan behavior. 

The distribution of Increments at lags 0.6, 6, 66, and 666 meters 
given in the Figures 7.5 - 7.8 all show the turbulence behavior from 
the Figures 6 . 3 and 6 . 4 with small lags deviating more sharply from 
Gaussian behavior than large lags. 
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Figure 7.1 Experiment 1. For further details see the legend 
for figure 6.1 
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Figure 7.2 Experiment 1. For further details see the legend 
for figure 6.2 
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Figure 7.3 Experiment 2. For further details see the legend 
for figure 6.1 
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Figure 7.5 Experiment 1, For further details see the 
legend for figure 6.3 


probability density probability density 


I 


223 


"Q 



-? -5 -a -1 D 1 2 5 4 F -F -4 -5 -2 -1 Q 12 5 4 5 

INCREMENTS L= 1DDD M INCREMENTS L= 1000 M 


Figure 7.6 Experiment 1. For further details see the 
legend for figure 6.4 
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Figure 7.7 Experiment 2. For further details see the 
legend for figure 6.3 
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Figure 7.8 Experiment 2. For further details see the 
legend for figure 6.4 
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The accumulated moments as given by the Figures 7,9 and 7.10 
show some differences from experiment to experiment to the actual 
turbulence. Figure 6.5, but no serious discrepancies. This is also 
the case for the exceddance plot. 

The auto~ and cross-correlations and the structure function. 

Figures 7.11 and 7.12, do not differ much between the experiments. 

From experiment to the actual turbulence. Figure 6.6, there is, 
however, one discrepancy; the experimental u and v fall off much too 
rapidly. 

The relation between the correlation function R(T) for a stationary 
time-series and the correlation function R^(T) calculated over the 
length T is given by 

(7.1) Kj,(T) = R(t) (1 - -^) 

The generated 5 sec pieces will on average have the same 
correlation function as the ensemble from which the eigenfunctions were 
obtained. This ensemble was obtained from the turbulence record, and 
although it is not justifiable to call this ensemble stationary, it 
seems that the transformation above can explain the main differences 
between Figure 6.6 and Figures 7.11 and 7.12, especially when it is 
remembered that the triangular window plots like an exponential function 
on log-linear axes. 






Figure 7.10 Experiment 2. For further details see the 
legend for figure 6.5 
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Figure 7 . 12 Experiment 2 . For further details see the 
legend for figure 6,6 
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The generating interval was 5 sec, and with a mean wind speed of 
-1 

6.6 m sec this corresponds to a length of 33 m, and so after a lag 
of 33 m the autocorrelation should drop to zero, which actually is 
the case. For the real turbulence, the zero value is reached around 
600 m. With this knowledge, we would expect the generated turbulence 
to look different from the turbulence recorded for the u and v components. 

The cross-correlations agree well for lags less than 30 m, and it 
is worth emphasizing the behavior of the important uw correlation. The 
remarks on the autocorrelation can also be applied to the structure 
function. 

The u, V, and w spectra are plotted on the Figures 7.13 - 7.15 to- 
gether with the turbulence spectra. The u and the v spectra fall off 
like the turbulence spectra with a -5/3 slope for frequencies larger 
than 1 Hz. Between 1 Hz and 0.03 Hz (scales 6.6 m to 330 m) , the 
epxerlmental spectra exceed the turbulence spectra, and from 0.03 Hz 
and less the opposite is true. We could have expected a pronounced 
peak in the experimental spectra at 0.2 Hz due to the generation of the 
turbulence in 5 sec pieces, but this is not the case; the effect, if 
any, has been spread over more than one decade of frequencies. 

The w-spectra follow the same pattern, although they seem to 
coincide much better except for high freqencles where the experimental 
spectra fall off too fast. This probably means that the eigenfunctions 
from 21 and up still have some significant high frequency features to 
add to w, but not to u and v. 

The implication of the transformation 7.1 would be a folding of the 

2 

spectra with a (sin oj/w) function, which except for some end effects 
will tend to preserve power law behavior. 
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The agreement between the uw-cospectra in Figure 7.16 and those 
in Figure 6.7 are very good. 

The whole turbulence time history, block averaged from 10 Hz to 
1 Hz, is displayed in Figure 6.8. The two experimental turbulence 
records were normalized to have the same means, zero, and variances 
as the turbulence record at 10 Hz and then block averaged to 1 Hz. 

As is evident from the former discussion, the differences betweeh 
the two experiments are small, and we will therefore further on 
confine our attention to experiment 1. Figure 7,17 shows 1000 sec 
of the experimental turbulence and as could be expected from the 
discussion of the correlation functions and the spectra, the u and v 
components have too much variation in the mid-frequency range, but 
apparently none at low frequencies. 

Figures 6.9 and 7.18 display 100 sec of the turbulence and the 
experiment respectively, all data at 10 Hz and normalized with mean 
zero and variance 1. As before, the w-components agree very well and 
so do the u and v components at Intervals less than 10 seconds. The 
effect of patching the 5 sec Intervals together is seen to cause some 
excessive large jumps, a deficiency which has to be corrected, possibly 
by establishing a patching procedure. 

7.2 Spectral Shaping 

In the final development of the model, some spectral shaping 
seems unavoidable, first because of the problems of generating energy 
at low wavenumbers, and secondly, because in practical applications 
care must be exercised so that the energy is maximized at the appropriate 


wavenumbers. 









d 



Figure 7.17 Experiment 1. Time history consisting of 
1000 second averages, normalized to the 
same variance as the turbulence record. The 
numbers on the horizontal and vertical axis 
denote seconds and meters per second 
respectively. 
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Several methods exist, probably the most well-'^.known being to 
let the spectra match the behavior of the von Karman or the Dryden 
spectrum (Smith, 1971, Fichtl, 1973). Another approach is to use 
the seml-emplrlcal spectral formulaes obtained by the mlcrmeteoro logical 
researchers (Busch, 1973). 

The fair degree of coincidence of the measured spectra aiid the 
experimental spectra ledds us to the conclusion that in order to shape 
the low wavenumber end of the experimental spectral we could as well 
shape them over the whole wavenumber range using the measured spectra. 
This was accomplished in the following way; 

Calculate the discrete Fourier transform for one component at a 
time and let the series of Fourier coefficients be given by 


%*^o»^l’^l 


the measured turbulence 


c d^ , c, , d- , . . . . c„ d„ the generated turbulence 
o, o X 1 N> N 


®o’^o’®l’^l’ 


the generated and shaped turbulence 


where 


(7.2) 


(ei,fi) = [c^ 


f 2 ^ . 2 . 1/2 

^ i ^i ^ 


(Ci + dj ) 


24/2 • 2 . ■ 24/2 ^ 


(c^ + ) 


which gives 
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(7.3) 

and 


f . d. 

1 X 



so that we have changed the spectral shape of the generated turbulence 
to that of the measured turbulence, but we have preserved the phase 
angles from the generated series. The series (e^,f^) is then back 
transformed to obtain the series in the time domain. 

The effects of the shaping are shown for experiment 1 in 
Figures 7.19 - 7.26, The autospectra and the autocorrelation are not 
shown as they obviously have to be the same for the measured record 
(the autocorrelation function approximately) . The agreement between 
the various statistical and sequential characteristics for the real 
turbulence and the simulated turbulence are evident. 

One desirable effect of the shaping is that it seemingly smooths 
out the effect of patching the generated 5 seconds together, and thereby 
relaxing the requirement for a patching procedure. 

7. 3 Summary Conclusions and Recommendations 

In this study several tasks were undertaken all with the primary 
goal of producing an operational turbulence simulation model. 

In Section 2, earlier attempts to use the Proper Orthogonal 
Decomposition in turbulence modeling were reviewed and the conclusions 
summarized in Table 1. The experience gained through these studies 
suggested the development of the present model. 
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3 figure 7.19 Experiment 1 after spectral shaping. For 
further details see the legend for 
figure 6.1 
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Figure 7.21 Experiment 1 after spectral shaping. For 
further details see the legend for 
figure 6.3 
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Figure 7.22 Experiment 1 after spectral shaping. For 
further details see the legend for 
figure 6.4 
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Figure 7.24 Experiment 1 after spectral shaping. Top: 

Co-spectra of two samples of the standardized 
u-and w-components. Each sample consists of 
8192 points (819.2 sec). Bottom: cross 

correlation for the standardized data plotted 
as functions of the lag on a logarithmic scale. 
The plotted 1 refers to the u- and v- 
correlation, the 2 to u-w, the 3 to v-w. The 
Integers on the horizontal axis denote 
power of 10 
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Because the Proper Orthogonal Decomposition method eventually 
has to extract information from one or more actual records, it was 
argued in Section 3 that under some justifiable assumptions of 
self similarity and spectral power- laws, this does not limit the 
generated turbulence to fixed length and time-scales. 

The theory behind the Proper Orthogonal Decomposition and its 
application in practice was briefly outlined in Section 4. The 
argument for applying this theory in the model went as follows. From 
earlier studies we know that it is essential to model the "surprise" 
in turbulence. The proper Orthogonal Decomposition provides us with 
a method to represent this phenomenon by a set of orthogonal functions 
in which the first function has the closest resemblence in a least- 
square sense to large gusts. The motion between these large gusts can 
also be represented by a set of orthogonal functions, which if provided 
by the Proper Orthogonal Decomposition, give a unique optimal ex- 
pansion of a stochastic process. Making the simplifying assumption 
that the motion in the gust intervals consists of a gust structure plus 
an orthogonal stochastic process allowed us to concentrate the 
analysis to certain selected gust intervals. 

The construction of the model was described as a two part process, 
first the analysis scheme to obtain the orthogonal functions, and 
second, the generating scheme which uses these functions. The two 
schemes were summarized at the end of Section 4. 

The results of the analysis scheme were discussed in Section 5 and an 
investigation of the eigenfunctions revealed the first few to be a 
mixture of gust structure and non-gust-motion with a decreasing emphasis 
on the gust as the order of the eigenfunctions Increases. The dis- 
tribution of the expansion-coefficients were estimated to be Gaussian. 
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The actual turbulence record was subjected to a fairly detailed 
analysis, reported in Section 6, in order to provide a set of statistics 
for comparison with those of the generated turbulence. The calculated 
statistics were compared with those obtained in other studies to ensure 
the representativeness of the record. The statistics were chosen with 
emphasis on the importance of the non-Gauss lan behavior of turbulence 
processes as well as the importance of the second order sequential 
statistics. 

Section 7 describes two experiments that were performed with 
the model and compared the generated turbulence with the actual 
turbulence. The analysis was performed twice before and after spectral 
shaping. The main conclusions before the shaping were: 

1. The two experiments differed no more than actual turbulence 
records. 

2. The non-Gaussian behavior was well-modeled. 

3. The sequentially dependent statistics and characteristics 
were modeled well on time scales less than the length of the eigen- 
functions . 

4. The time history showed some excessive jumps due to patching 
generated intervals together. 

The spectral shaping was accomplished by transforming the gen- 
erated turbulence into series with spectra equal to those of the actual 
turbulence, and preserving the phase angles during the process. After 
the spectral shaping, the experimental turbulence appeared to be as 
close to real turbulence as any practical model might produce. 

However, the model is not perfect and for future work we would 
recommend that: 
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1. More, turbulence records be used in the analysis. 

2. The length of the ensemble functions selected to determine 
the eigenfunctions to be chosen with more considerations to the 
integral scales of the process, 

3. Several choices of the feature to be tested, for example 
certain events in the uw and uvw correlations, 

4. The probability density function and the time interval* be- 
tween the events should be estimated properly from the turbulence 
records , 

5. The effects of preparing the selected ensemble in different 
ways should be investigated, 

6. The difficult choice of the transformation of the density- 
functions for the first eigenfunctions could be eased by extending 

the analysis so that two sets of eigenfunctions are calculated, one set 
for the integrals with special features and one set for all other 
Intervals. Both sets could then be used in the generating scheme, or 
they could be mixed, depending on an investigation of the similarities 
between the two sets. This procedure would provide a test of whether 
a characteristic feature exists in the chosen ensemble, because if the 
two sets of eigenfunctions show close resemblance, this is clearly not 
the case. 

7. Alternative spectral shapings should be Investigated. 

The actual computer programs that generate the turbulence are 

rather simple and not very time consuming, and so a lot of experiments 
can easily be performed, in order to establish how the properties of 
the generated turbulence vary with changes in the model. However, it is 
a laborious task to make a thorough study of each experiment as can be 



252 


judged from the number of figures in this report. A solution might 
be to investigate just a few experiments with extreme variations in 
the parameters, and then to test other experiments by applying the 
generated turbulence in practice and judging the outcome. 
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8.0 SOME CONSIDERATIONS ON USE OF THE KARHUNEN-LOEVE 
EXPANSION IN DATA ANALYSIS 

Erik L. Petersen 

The Proper Orthogonal Decomposition Theorem has been used by an 
increasing number of researchers in a variety of fields to Investigate 
time or space series for quasi-deterministic structures. Among them 
are Lorenz (1956), Lumley (1965), Dutton, et al, (1968), (1969), (1971), 
Holstr^m (1970), Jasperson (1971), Busch and Petersen (1971), just to 
mention a few. 

Let us here be concerned with whether quasi-deterministic behavior 
appears in the functions in (f(t)} where {f(t)} is an ensemble of second 
order real valued random functions of the parameter -oo < t < °° and where 

E (f(t)} = 0 

(8.1) E {f^(t)> = 1 

E {f(s)f(t)} = R(s,t) 


Let {f^} be subensembles formed from {f} by assembling the sequential 
values over Intervals of length T of some functions f£,{f}. For example, 
if 


(8.2) H(x) 



then H(x+a)f(x) would be in one of the ensembles {f^} for every a. For 
ease of computation, each function in 'Cf.j,} is redefined over the domain 
[0, T]. 
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If we suppose a quasi-determlnistlc component appears in {f}, 
then the question is how to select a criterion that will give a 
function 4>(t), 0 £ t ^ T, that resembles the deterministic part of 
the functions in an optlmtim manner. 

Several measures of resemblance are possible; the one we will 
choose is quadratic and has been discussed previously (Lumley, 

1965, Dutton, 1969). 

Let {f^} be one of the subensembles of {f}. Then we define 

fT 

[ ^ f(t) <Mt) dt]^ 

(8.3) A = Ef , 

The function (f) we are seeking maximizes A over the collection of all 
subensembles {f^} for various values of T. Adopting other criteria 
in order to determine the quasl-deterministic behavior given by 4> 
would in general lead to other approaches, but it is the criterion 
above that brings the Karhunen-Loeve expansion into the analysis. This 
expansion known to have some very general properties, and this in turn 
justifies the criterion. 

The question whether there exists a unique solution to (8.4) such 
that it is possible to find a subensemble {f^} determining a function cf)(t) 
which gives an absolute maximum of A will be considered in the last 
part of this section. 

If a maximum is found, let {f^}* be the associated ensemble and let 
{f^} be constructed as 

{f.^r - 
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which is to say that are subensembles formed from {f} by 

assembling sequential values over intervals of length T which have not 

■ r n * 

been used in the construction of if^j . 

We will now make the assumption that almost all the information 
contained in {f} can be estimated from <J)(t) and {f^} . A necessary 
but not a sufficient condition for this assumption to be true is 
that R(s,t) does not differ significantly from zero outside |s-tj > T. 
Let us further assume that we need a representation of the sequential 
characteristics of {f^} which is as economical as possible. Such a 
representation may be found by expanding {f^} in complete orthogonal 
systems, if it is possible to find an expansion that gives a good 
approximation to the ensemble. functions by an economically small number 
of terms. 

The optimal expansion of a random function, the Karhunen-Loeve 
expansion, is suitable for this purpose. The expansion is optimal 
in the sense that the series truncated at any point minimizes the 
integrated mean square deviation between the actual and the approximated 
random functions. Any other expansion using the same number of terms 
cannot have an Integrated mean square deviation which is less. This 
is to say that minimizing the error e(N) 

f N 9 

E{e(N)} = E{ |f(t) - Z a ij; (t) | dt} 

J ^ n n ' 

^ n=l 

leads to an expansion (see e.g. Dutton, 1969) 

00 

f(t) = Z a il^.(t) 
k=l ^ ^ 


(8.4) 
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where the functions used in this expansion are the eigenfunctions 

of the Fredholm integral equation 

(8.5) j R(s,t) dt = 

o 

and 

(8.6) J i|)^(t) iJ)jj^(t)dt = 6 ^^ 

(8.7) E {a a } = A 6 

n m n nm 

(8.8) a = f(t) ^„(t) dt 

n J n 

where the X^*s have been arranged in a non’- increasing sequence. 

This follows from the Proper Orthogonal Decomposition Theorem, 

Loeve (1955) , which states that a mean’-square continuous random 
function f(t) defined on a closed interval 0 ^ t _< T, has the 
decomposition (8.4) with the properties given by (8.6), (8.7), and 

(8.8) if and only if X^ are the eigenvalues and 'l^^(t) the orthonormal 
eigenfunctions belonging to the correlation function R(s,t), and thus 
are solutions to (8.5). The theory is based on Mercer's Theorem, which 
states that a non-negative definite function, R(s,t) continuous over 
the closed interval 0 ^ s,t ^ T has the expansion 

(8.9) R(s,t) = EX ^„(s) ^„*(t) 

n n n 
n 

where X^ and i|^jj(t) are the solutions of (8,5) ^nd the asterisk denotes 


a complex conjugate. 



Solving the equations (8.5) - (8,8) gives us the Karhunen’-Loeve 
expansion, sometimes referred to as a generalized spectral represen- 
tation because of the lack of correlation between the coefficients 
of expansion. 

From the equations (8.7) and (8.8) we get 



which shows that the eigenvalues reveal the fraction of the total 
variance which Is explained by the associated eigenfunction. 

Our problem now is to find {f,j,} , and thereby {f^} and (j), by 
solving the variational problem as given by (8.3). Unfortunately, we 
are not able to do this In full generality. We cannot analytically 

r 1* 

find the subensemble If^j that maximizes A over all subensembles 
{f^}; however. If by some other methods we can establish {f^} we can 
solve (8.3) for the function (f) which maximizes X over the ensemble 
If^J- . An approximate solution could be obtained by calculating X and 
plotting It for various choices of {f^}, but the work Involved Is 
staggering. 

Let us then assume that we are able to select a subensemble 
{f.j,} subjectively which Is not far from the optimizing subensemble 
If one exists. 

Applying the techniques of the calculus of varltlons to (8.3) to 
find the maximizing function, (j)(t), leads to the Integral equation; 
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and we observe that the function most like each of the functions in the 
ensemble {f^} is determined by the Karhunen'-Loeve expansion of {f^}, and 
we have <f)(t) = (})^(t). The expansion can be interpreted in the following 
way: 4)^^ is the single function that explains the most variance in the 

ensemble, but it does not explain all of the variance, so we form a new 
ensemble of functions {f^ - a^ a^id find the one function most like 
the residual function; the answer will be and so we consider a new 
ensemble {f^ - a^ - a^ 

With the assumptions made through this section we are now able to 
represent the sequential characteristics of {f} by the two orthonormal 
systems, and together with the sampling properties of the cor- 
responding expansion coefficients. 

The two systems can be reduced to one system if further assumptions 
are made: if the quasideterministlc structure we are seeking occurs 

over the length T, and if this structure is orthogonal to the rest 
of the process taking place over the length T, and if this structure is 
given entirely by the first eigenfunction then the ensemble 
{f^ - a^ 4)^} will have the same properties as {f,j,} hence 4*^ will be 
the same as ^2 assumption was made in the 

construction of the model discussed in the preceding part 
of this report. 

But does there exist a unique solution to (8.3) such that it is 
possible to find a subensemble {f^} determining a function 4>(t) which 
gives an absolute maximum of A? An example will show that such an 
ensemble does not exist in general. 
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8.1 The Karhunen-I^eye Expansioii of a First Order Autoregressive Series 
For the example we will use a stationary time series f (t) with the 

I tl^S I 

correlation function R(s~t) = e ' , a > 0. This Is the correlation 

function for a first order autoregressive series. 

f(t) = R(t» 1) f(t-l) + e(t), t = p, 1, 2, ... 

where e(t) Is normally distributed with E{e(t)} = 0 and 

2 

E{e(t) e(s)} =6 o , This form of a correlation function Is 

S ^ u c 

quite often a good approximation for geophysical time series despite 
Its lack of microscale, l.e., R(t) Is not differentiable for T = 0. 

A very obvious reason for choosing R(t) Is that It Is then possible 
to find an analytic solution to the Fredholm Integral equation (8.5). 

Let us first prove a theorem that often can be useful. 

Theorem: The Karhunen-Loeve expansion of a stationary second’-order 

process yields eigenfunctions that are either odd or even. 

Proof: From Mercer’s Theorem, (8.9) we have 

(8.10) R(s-t) = Z A ^ ^ s,t ^ T 

n 

which upon replacing s with T-s and t with T-^t gives 

(8.11) R(-s+t) = Z A^ 6 (T-s) 6 (T-t) 

n n n 

n 

But R Is an even function so that R(-s+t) = R(s-t) , and thus combining 


(8.10) and (8.11) gives 
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R(s-t) = S X (j) (s) 4>_(t) = E X <(> (T-s) A (T-t) 
n n II n n n 

n n 

From the uniqueness of the Karhunen-Lbeve expansion It then follows 
Immediately that for every 0 ^ t £ T 

0„(t) = + b (T-^t) 
n — n 

so that the eigenfunctions are either odd or even with respect to the 
point T/2. This concludes the proof. 

Although the converse Is not true because a set of odd or even 
eigenfunctions only ensure that we have a correlation matrix symmetric 
with respect to the two main diagonals, one could use the theorem with 
care to estimate the non-statlonarlty of a certain ensemble. If, for 
example, the eigenfunctions appear to be even or odd from the m^^ 
eigenfunction on, arid If the first m-1 eigenfunctions explain P% of the 
variance, then a reasonable hypothesis would be that approximately P% 
of the variance in the ensemble could be due to Instatlonarltles. 

Proceeding with the example, let us consider an ensemble over the 
closed Interval 0 t ^ T with the following properties 

E {i j f (t) dt} - 0 

o 

7 

E J f^(t) dt} = 1 
o 

E {i j f(t) f(t-T) dt} = 
o 
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The eigenfunctions of the correlation function are determined by 


( T 

^-a|t-s| ^ (s)ds ■ X 4) (t) 
^n n ^n 


or, dropping the subscript n 


r 


(8.12) X4>(t) * I e”^^^ 4>(s)d8 + 1 e (j>(s)ds 


-a(s-t) 


Differentiating twice, using the rule for differentiating under the 
Integral sign, we have first 


X4»'(t) 


-•r 


e~^(t ({)(8)ds + oi 


rT 


-a(s-t) .. 

e <p(8)ds 


and then 


X(j) 


’*(t) = a^[f^ e‘""^^"®^(|)(s)ds + r e"^^®“^^(t)(s)ds]-2a(|>(t) 
Jo t 


so that 


4>"(t) + ~ ^ <C(t) “ 0 


As a first case, we consider 


m r. g...A< 0 -> X > -^ 

X a 


The general solution to the differential equation is 
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♦ (t) - + C2e"®* 


but It is easily seen that no colce of C^, C 2 can make <{>(t) either 
odd or even over the Inter-^^al (0, T) hence, no solutions exist for 
\ > 2/a. 

As a second case, we have A < 2/a. The general solution Is now 


4>(t) * cos(a)t) + sin(wt), o) = 



or 


(|)(t) =* A sln(a>t + 0j|^) 
which we can write as 


(f)(t) »» A sin (u)(t " + 0) 


To determine A, O), 0 and A we have the Integral equation, (8.12), 
and the condition of orthogonality (8.6). After some algebra we find 

(8.13) i|)^(t) » /y sin (U)^(t - |) + 22 :) 

n 

and the eigenvalue spectrum is given by 


(8.14) A = 
n 


2a 


2 ^ 2 
a + 0 ) 

n 
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where the are solutions to the equation 


2a oj 


(8.15) tan(o) T) = 
n 


n 


2 2 
a -^ 0 ) 
n 


where 


(8.16) “n ^ “ = 


The normalization of the ensemble gives 


(8.17) Z A = 1 

i n 
n 


The integral scale for the process f(t) is defined as 


I = 


T 1 

e ’ ' dT = — 

a 


A variable C is introduced by 


(8.18) C = aT (= Y> 


which is seen to be the ratio between the chosen length of the eigen^ 
functions and the integral scale. 

Introducing C into (8.14) and (8,15) yields; 


(8.19) 


n 


2C 


T 2 2 

^ C + (w T) 
n 


I 
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2C (o) T) 

(8.20) tan (w T) = - 

" C - (a> T) 

n 


We can now numerically solve (8.17) ~ (8,20) to give n - 

1, 2, as a function of C. The result is shown on Figure 8.1 

on a logarithmic scale. 

As expected, it is not possible to find, in a statistical sense, 
r 1 * 

an ensemble if^j that gives an absolute maximum for because 
X^/T -> 1 for T, C 0. 

It is thus seen that applying the Karhunen-Loeve expansion to a 
certain selected ensemble and finding that the first eigenfunction 
accounts for a tremendous amount of the variance does not automatically 
ensure that a characteristic structure is revealed by the first eigen- 
function. 

Before we conclude this example, we will show an interesting 
relationship, in this case, between the Fourier spectrum and the 
eigenvalue spectrum. 

The normalized Fourier spectrum of f(t) is given by 


(8.21) S(03) = ^ 


-ax -icoT , 2 

e e dx = — 


a 


7T 2 ^ 2 

a + w 


and the normalized eigenvalue spectrum by 


X 

n 


n n 


g 

T 2 ^ 2 

g + 0) . 


(n-l)7T 

T 


CO 

n 



and then all the eigenvalues are, except for a constant factor, lying 
on the curve given by the Fourier spectrum. 



RATIO OF VARIANCE 



n . 01 


0. 10 I. 0 

TIMESCALE / INTECRALSCALE 


Figure 8.1 Curves showing for each of the first six 
eigenfunctions of a first order auto- 
regressive series the amount of normalized 
variance explained as function of the ratio 
of the length of the eigenfunctions and 
the integral scale of the series. The 
curve for eigenfunction n+1 is always below 
the curve for eigenfunction n 
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The mean value theorem gives 


IT 


(n-l)Y 


a 


, 7T 2 a _ 2 

S (0)) d(i) = =: =• TT — 7z 

Ttt2^ 2 T2^ 

ct + u)£ a ^1 


_ (n-l)TT < < n# 

2* T A T 


where O) for T ^ 

Jc n 

Two approximate 'formulas: 


n ^ 2 

— - ^rc tan 


Ctt 


+ n (n-1) TT^ 


m A 

■ li . ^ ^ Mill 

^ t" " T c“ 

n~l 


for Wfl 03 
Si n 


8.2 The Karhunen-Loeve Expansion of a Bandllmlted White Noise Process 

A slightly different approach for solving the Integral equation (8.5) 
with an exponential kernel can be found In Davenport and Ross (1958) and 
Pugachev (1965) and (1959). In Sleplan, Poliak, and Landau (1961) a 
specially Interesting case Is analyzed where the tlmeserles Is band- 
llmlted white noise: 


S(03) * 


1/2 Q, |o)| < 

0 , |o3| > 


“ I In 

-fi 


sin fit 
At 


The solutions to 
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T/2 


sin Q(t-s) 

n(t-s) 


4» (s) 

n 


ds 


X <^ (t) 
n 


-T/2 

are knovn as prolate spheroidal wave functions ip 

c - fiT 


n 


The eigenfunctions and the eigenvalues depend only on c or say 


- [ R(t) dT = 
•' o 


sin At _ 1 2L 
nr 2 


c 


^^T 


n ^ 
2 1 


They depend as in the previous example only on the ratio between the 
length of the ensemble function and the Integral scale of the process. 


8,3 Conclusion 

The Karhunen-tioeve expansion (or the Proper Orthogonal Decomposition 
or the Generalized Spectral Representation) has often proved to be a 
powerful tool in mathematical-statistical analysis of random processes 
especially where advantages can be taken of the lack of correlation 
between the expansion coefficients. However, the method Is not without 
disadvantages when coiiq>ared with the usual Fourier analysis. One is 
that the eigenfunctions have to be tabulated or plotted while sines and 
cosines are well-known functions. Another Is that a physical Interpre- 
tation of the eigenvalue spectrum Is often difficult, if possible, at all. 
In contrast to Fourier spectra where the amount of explained variance 
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Is plotted versus frequency or wavenumber and hence makes direct 
reference to time or length scales In the process being studied « 

Fourier analysis usually requires the process to be stationary 
and ergodic whereas the Proper Orthogonal Decomposition can be 
applied to non-s tat ionary ensembles. But in order to use the 
latter method, it is necessary to establish an ensemble, and in the 
process of doing so it is of the utmost importance to be sure that all 
the ensemble functions included are equal members of the ensemble. 

It is obvious that in order for the method to reveal a characteristic 
structure in the ensemble, the ensemble functions have to be properly 
aligned. If for example the structure is a sinusoid, there must be no 
phaseshift from one ensemble function to another. 

In this section we have outlined the theoretical basis for the 
model which was developed in the preceding sections. We believe that 
this basis can support further Investigations of the potential of the 
Proper Orthogonal Decomposition in studying and modeling physical 
processes. We have also pointed to the necessity of investigating the 
statistical structure of the series being studied before it is decided 
how to create the ensemble functions the method actually required as 
input . 
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